# a vector of all the packages needed in the project
packages_required_in_project <- c("RMark",
"tidyverse",
"readxl",
"BaSTA",
"pbapply",
"RColorBrewer",
"grid",
"Rmisc",
"gss",
"arm",
"partR2",
"parameters",
"standardize",
"colorBlindness",
"ggthemes",
"patchwork",
"gt",
"rptR",
"tidybayes",
"broom.mixed",
"effects",
"patchwork",
"devtools",
"unmarked",
"R2ucare",
"marked",
"merTools",
"bootpredictlme4",
"extrafont",
"survminer",
"bdsmatrix",
"coxme",
"epitools",
"survival",
"magrittr",
"ggpubr",
"reshape2",
"sp",
"adehabitatLT",
"reshape",
"coefplot",
"mapview",
"lubridate",
"effects",
"merDeriv",
"remotes")
# of the required packages, check if some need to be installed
new.packages <-
packages_required_in_project[!(packages_required_in_project %in%
installed.packages()[,"Package"])]
# install all packages that are not locally available
if(length(new.packages)) install.packages(new.packages)
# install the bootpredictlme4 package directly from GitHub
remotes::install_github("RemkoDuursma/bootpredictlme4")
# load all the packages into the current R session
lapply(packages_required_in_project, require, character.only = TRUE)Vignette of analysis: Early-life demographic processes do not drive adult sex ratio biases and mating systems in sympatric coucal species
Supplementary Material A
In this document we provide all the necessary code for reproducing the analyses presented in our manuscript. To access the dataset and Rmarkdown file, please download this GitHub repository. Simply follow the link and click on Download ZIP on the right-hand side of the page. An explanation of the files in the repository can be found in the Readme file. Please don’t hesitate to contact Luke at luke.eberhart[at]bi.mpg.de if you have any questions.
Prerequisites
GitHub repository
- For running the complete code you need a
filessubfolder containing the raw data downloaded fromData_filesfolder provided in the GitHub repository.
R packages
- The following packages are needed for analysis and can be easily installed from CRAN using the following code:
Plotting themes
- The following plotting themes, colors, and typefaces are used throughout the project:
# Find fonts from computer that you want. Use regular expressions to do this
# For example, load all fonts that are 'candara' or 'Candara'
extrafont::font_import(pattern = "[V/v]erdana", prompt = FALSE)
# check which fonts were loaded
extrafont::fonts()
extrafont::fonttable()
extrafont::loadfonts() # load these into R
# define the plotting theme to be used in subsequent ggplots
luke_theme <-
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.title = element_text(size = 16),
legend.text = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 10),
axis.title.y = element_text(size = 12),
axis.text.y = element_text(size = 10),
strip.text = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size = 0.5, colour = "grey40"),
axis.ticks.length = unit(0.2, "cm"),
panel.border = element_rect(linetype = "solid", colour = "grey"),
legend.position = c(0.1, 0.9)
)
# set plotting color palettes
sex_pal2 <-
c(pull(ggthemes_data$wsj$palettes$colors6[3,2]),
pull(ggthemes_data$wsj$palettes$colors6[2,2]))
sex_pal3 <-
c(pull(ggthemes_data$wsj$palettes$colors6[3,2]),
pull(ggthemes_data$wsj$palettes$colors6[3,2]),
pull(ggthemes_data$wsj$palettes$colors6[2,2]),
pull(ggthemes_data$wsj$palettes$colors6[2,2]))
plot_palette_sex <- RColorBrewer::brewer.pal(8, "Dark2")[c(2,1)]
plot_palette_species <- RColorBrewer::brewer.pal(8, "Dark2")[c(3,7)]
# specify the facet labels for each species and sex
species_names <- c(
'BC' = "Black coucal",
'WBC' = "White-browed coucal")
sex_names <- c(
'female' = "Females",
'male' = "Males")
# color of mean estimate point in forest plots
col_all <- "#2E3440"Analytical Functions
The following custom functions are used to bootstrap the various datasets and run the stochastic demographic model. To reproduce the results presented in our manuscript, these functions need to be stored in the environment of your R Session.
adult_CJS_bootstrap_data()
This function randomly samples the adult mark-recapture dataset to produce a random one the same size as the original, with replacement
arguments:
- adult: the adult mark-recapture dataset
- num_boot: a unique number for the pbapply simulation (don’t specify)
- species: species name (either “BC” or “WBC”)
- iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
adult_CJS_bootstrap_data <-
function(adult, num_boot, species, iter_add) {
# sample a new adult mark-recapture dataset the same size as the original,
# with replacement
adult_boot <-
adult %>%
sample_n(nrow(.), replace = TRUE)
# make a list of the new dataset, iteration, species
out <- list(adult_boot = adult_boot,
iter = num_boot + ((iter_add - 1) * niter),
species = species)
}bootstrap_adult_survival_CJS()
runs the CJS survival analyses using the the bootstrapped sample created from adult_CJS_bootstrap_data().
arguments:
- coucal_boot_list: the bootstrapped adult mark-recapture dataset produced from the adult_CJS_bootstrap_data() function above
- num_boot: a unique number for the pbapply simulation (don’t specify)
- first_year: first year of the study (i.e., the year of the first event in the encounter history)
- bootstrap_name: a unique string for naming of files in the output directory
- species: species name (either “BC” or “WBC”)
- iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
- prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
- out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
bootstrap_adult_survival_CJS <-
function(coucal_boot_list,
num_boot,
first_year,
bootstrap_name,
species,
iter_add,
prefix_number,
out_directory = "/output/bootstraps/raw/") {
# specify the bootstrapped data samples (from the previous function)
adult_ch <- coucal_boot_list[["adult_boot"]]
# clean up capture histories
adult_ch <-
adult_ch %>%
ungroup() %>%
as.data.frame() %>%
mutate(sex = as.factor(sex))
# process the adult data as a "CJS" analysis
coucal_adult.proc <- RMark::process.data(adult_ch,
model = "CJS",
groups = c("sex"),
begin.time = first_year)
# Create the design matrix from the processed mark-recapture datasets
coucal_adult.ddl <- RMark::make.design.data(coucal_adult.proc)
adult_survival_analysis_run =
function(proc_data, design_data){
# sex- and stage-specific survival:
Phi.sex = list(formula = ~ sex)
# Models exploring variation in encounter probability
# constant:
p.dot = list(formula = ~ 1)
# sex-dependent:
p.sex = list(formula = ~ sex)
# factorial variation across year:
p.Time = list(formula = ~ time)
# additive effects of sex and factorial year:
p.sex_Time = list(formula = ~ sex + time)
# create a list of candidate models for all the a models above that begin with
# either "Phi." or "p."
cml <- RMark::create.model.list("CJS")
# specify the data, design matrix, delete unneeded output files, and
# run the models in Program MARK
model.list <- RMark::mark.wrapper(cml, data = proc_data,
ddl = design_data, delete = TRUE,
wrap = FALSE, threads = 1, brief = TRUE,
silent = TRUE, output = FALSE, prefix = prefix_number)
# output the model list and sotre the results
return(model.list)
}
# Run the models on the bootstrapped data
adult_survival_analysis_out <-
adult_survival_analysis_run(proc_data = coucal_adult.proc,
design_data = coucal_adult.ddl)
extract_top_model_output <-
function(rmark_output, top_model = TRUE, mod_num){
# Find the model number for the first ranked model of the AIC table
if(top_model == TRUE){
mod_num <-
as.numeric(rownames(rmark_output$model.table[1,]))
}
else{
mod_num <- mod_num
}
# extract and wrangle reals from model output
reals <-
rmark_output[[mod_num]]$results$real %>%
bind_cols(data.frame(str_split_fixed(rownames(.), " ", n = 5)), .) %>%
filter(X1 == "Phi") %>%
mutate(sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F", "Female", "Male"))) %>%
dplyr::select(sex, estimate) %>%
mutate(iter = num_boot + ((iter_add - 1) * niter),
species = species)
# extract and wrangle the beta estimates from the model output
betas <-
rmark_output[[mod_num]]$results$beta %>%
mutate(statistic = rownames(.),
iter = num_boot + ((iter_add - 1) * niter),
species = species) %>%
mutate(parameter = ifelse(grepl(x = statistic, pattern = "S:"), "S",
ifelse(grepl(x = statistic, pattern = "p:"), "p",
ifelse(grepl(x = statistic, pattern = "r:"), "r",
ifelse(grepl(x = statistic, pattern = "F:"), "F",
ifelse(grepl(x = statistic, pattern = "Phi:"), "Phi", "XXX"))))),
variable = ifelse(grepl(x = statistic, pattern = "Intercept"), "Intercept",
ifelse(grepl(x = statistic, pattern = "sexM:Cubic"), "sexM:Cubic",
ifelse(grepl(x = statistic, pattern = "sexM:Quadratic"), "sexM:Quadratic",
ifelse(grepl(x = statistic, pattern = "sexM:Time"), "sexM:Time",
ifelse(grepl(x = statistic, pattern = "sexM"), "sexM",
ifelse(grepl(x = statistic, pattern = "Time"), "Time",
ifelse(grepl(x = statistic, pattern = "Cubic"), "Cubic",
ifelse(grepl(x = statistic, pattern = "Quadratic"), "Quadratic","XXX"))))))))) %>%
dplyr::select(-statistic)
# consolidate the AIC model selection results
AIC_table <-
rmark_output$model.table %>%
mutate(iter = num_boot + ((iter_add - 1) * niter),
species = species,
model_no_orig = as.numeric(rownames(.))) %>%
mutate(model_no_rank = as.numeric(rownames(.)))
# consolidate the all output into a list
survival_model_output_list <-
list(reals = reals,
betas = betas,
AIC_table = AIC_table)
survival_model_output_list
}
# extract and format survival rates from model output
adult_survival_model_output_list <-
extract_top_model_output(rmark_output = adult_survival_analysis_out)
# make a list of all the results from this iteration
bootstrap_results_list <-
list(adult_survival_model_output = adult_survival_model_output_list,
bootstrapped_data = coucal_boot_list)
# extract file directory
file_directory <- paste0(getwd(), "/", out_directory, "/", bootstrap_name, "_", num_boot, ".Rds")
# save the results of the iteration
save(bootstrap_results_list, file = file_directory)
bootstrap_results_list
}run_bootstrap_adult_survival_CJS
run the bootstrap and the CJS analysis (passed through the pbsapply function), see arguments defined abov
arguments:
- adult: the adult mark-recapture dataset
- num_boot: a unique number for the pbapply simulation (don’t specify)
- species: species name (either “BC” or “WBC”)
- iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
- first_year: first year of the study (i.e., the year of the first event in the encounter history)
- bootstrap_name: a unique string for naming of files in the output directory
- prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
- out_directory: location in the project directory where the output files should be stored (e.g., “tmp”)
run_bootstrap_adult_survival_CJS <- function(adult, num_boot, species, iter_add,
first_year, bootstrap_name, prefix_number,
out_directory)
{
# run the sampling function and specify the datasets
bootstrap_data_list <-
adult_CJS_bootstrap_data(adult = adult,
num_boot = num_boot,
species = species,
iter_add = iter_add)
# run the survival analysis and ASR deduction on the sampled data
result <- bootstrap_adult_survival_CJS(coucal_boot_list = bootstrap_data_list,
num_boot = num_boot,
first_year = first_year,
bootstrap_name = bootstrap_name,
species = species,
iter_add = iter_add,
prefix_number = prefix_number,
out_directory = out_directory)
result
}surv_boot_out_wrangle()
loads all the output rds files from the parallel adult CJS bootstrapping procedure and consolidates them into a single object for further analysis
arguments:
- species: species name (either “BC” or “WBC”)
- niter: number of iterations specified in the stochastic simulation (used to predict how many files to merge from the output directory)
- out_directory: location in the project directory where the output files from the simulation were stored (e.g., “tmp”)
# boot_out_wrangle()
surv_boot_out_wrangle <- function(species, niter, output_directory = "output/bootstraps/"){
# specify directory string
boot_out_path <- paste0(output_directory,
species,"_adult_survival_CJS_bootstrap_result.rds")
# load the bootstrap output file
survival_bootstrap_result <- readRDS(boot_out_path)
# bind all the adult rates output together
adult_reals_survival_rates_boot <-
rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$reals))) %>%
arrange(iter, sex)
# bind all the adult beta estimates output together
adult_betas_survival_rates_boot <-
rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$betas)))
# bind all the adult AIC table output together
adult_AIC_tables_boot <-
rbind(do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
function(x) survival_bootstrap_result["adult_survival_model_output", x]$adult_survival_model_output$AIC_table)))
# tidy all the bound data together into a mega list of results
boot_out_tidy <-
list(adult_reals_survival_rates_boot = adult_reals_survival_rates_boot,
adult_betas_survival_rates_boot = adult_betas_survival_rates_boot,
adult_AIC_tables_boot = adult_AIC_tables_boot)
boot_out_tidy
}approx_sd()
derive the approximate sd given the upper and lower 95% confidence intervals (see:Cochrane Chapter 7.7.3.2)
arguments:
- x1: lower bound of the confidence interval
- x2: upper bound of the confidence interval
# function to approximate the sd given the 95% CIs
approx_sd <- function(x1, x2){
(x2-x1) / (qnorm(0.975) - qnorm(0.025))
}matrix_ASR()
calculates the ASR of the population based on the two-sex two-stage projection matrix built by the coucal_matrix() function (below).
arguments:
- M: a two-sex x by x projection matrix produced by the coucal_matrix() function above
- n: x-lengthed vector representing starting stage distribution (the default is a vector with 10 individuals in each stage)
- h: harem size index (default is 1, monogamy; polyandry < 1, polygyny > 1)
- k: clutch size (default is 4 eggs, mean for both coucal species)
- iterations: number of time steps to project the demogrphic model (n.b., after 100 steps our model reaches the stable stage distribution)
- HSR: hatching sex ratio (expressed as the proportion of male chicks)
- num_boot: a unique number for the pbapply simulation (don’t specify)
- species: species name (either “BC” or “WBC”)
- iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk (if you run this on one session, then don’t specify)
matrix_ASR <-
function(M,
n = rep(10, nrow(M)),
h = 1,
k = 4,
iterations = 1000,
HSR = 0.5,
num_boot,
species,
iter_add){
# Number of stages in matrix
x <- length(n)
# Number of time steps to simulate
t <- iterations
# an empty t by x matrix to store the stage distributions
stage <- matrix(data = numeric(x * t), nrow = x, ncol = t)
# an empty t by x matrix to store the stage-specific frequencies
pop_dist <- matrix(data = numeric(x * t), nrow = x, ncol = t)
rownames(pop_dist) <- rownames(M)
colnames(pop_dist) <- 0:(t - 1)
# an empty t vector to store the population sizes
pop <- numeric(t)
# for loop that goes through each of t time steps
for (i in 1:t) {
# stage distribution at time t
stage[,i] <- n
# population size at time t
pop[i] <- sum(n)
# number of male adults at time t = (number alive at t-1) * (survival rate)
M2 <- stage[4, i] * M["M_Adult", "M_Adult"]
# number of female adults at time t
F2 <- stage[2, i] * M["F_Adult", "F_Adult"]
# Female freq-dep fecundity of Female chicks
M["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - HSR))
# Female freq-dep fecundity of Male chicks
M["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)
# Male freq-dep fecundity of Female chicks
M["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - HSR))
# Male freq-dep fecundity of Male chicks
M["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)
# define the new n (i.e., new stage distribution at time t)
n <- M %*% n
# define rownames of stage matrix
rownames(stage) <- rownames(M)
# define colnames of stage matrix
colnames(stage) <- 0:(t - 1)
# calculate the proportional stable stage distribution
stage <- apply(stage, 2, function(x) x/sum(x))
# define stable stage as the last stage
stable.stage <- stage[, t]
pop_dist[, i] <- n
}
# calc ASR as the proportion of the adult stable stage class that is male
ASR <- stable.stage[4]/(stable.stage[2] + stable.stage[4])
SSR <- stable.stage[3]/(stable.stage[1] + stable.stage[3])
# make a list of results
pop.proj <- list(ASR = ASR,
SSR = SSR,
lambda = pop[t]/pop[t - 1],
pop_size = pop,
stage_size = pop_dist,
stable.stage = stable.stage,
stage.vectors = stage,
SSD_M2 = stable.stage[4],
SSD_F2 = stable.stage[2],
iter = num_boot + ((iter_add - 1) * niter),
species = species)
# print the list as output to the function
pop.proj
}coucal_matrix()
builds the two-sex Lefkovitch matrix using the vital rates specified in the demographic_rates object. This function also give the option to specify a one-sex matrix for analytical comparisons between 2- and 1-sex matrix models
arguments:
- demographic_rates_: a list of all the vital rates of an iteration of the stochastic simulation run within the run_bootstrap_juv_hazd_ad_surv_ASR() function below
- two_sex: whether to make a one-sex or two-sex matrix of the rates provided (default is ‘two_sex = TRUE’)
coucal_matrix <-
function(demographic_rates_, two_sex = TRUE){
if(two_sex){
# Define coucal life-stages of the coucal matrix model
stages <- c("F_1st_year", "F_Adult", "M_1st_year", "M_Adult")
# Build the 4x4 matrix
result <-
matrix(c(
# top row of matrix
0, NA, 0, NA,
# second row of matrix
(demographic_rates_$Egg_survival *
demographic_rates_$F_Nestling_survival *
demographic_rates_$F_Groundling_survival *
demographic_rates_$F_Fledgling_survival),
demographic_rates_$F_Adult_survival,
0, 0,
# third row of matrix
0, NA, 0, NA,
# fourth row of matrix
0, 0,
(demographic_rates_$Egg_survival *
demographic_rates_$M_Nestling_survival *
demographic_rates_$M_Groundling_survival *
demographic_rates_$M_Fledgling_survival),
demographic_rates_$M_Adult_survival),
nrow = length(stages), byrow = TRUE,
dimnames = list(stages, stages))
}
else{
# Define coucal life-stages of the Ceuta snowy coucal matrix model
stages <- c("1st_year", "Adult")
# Build the 4x4 matrix
result <-
matrix(c(
# top row of matrix
0, NA,
# second row of matrix
(demographic_rates_$Egg_survival *
demographic_rates_$Nestling_survival *
demographic_rates_$Groundling_survival *
demographic_rates_$Fledgling_survival),
demographic_rates_$Adult_survival),
nrow = length(stages), byrow = TRUE,
dimnames = list(stages, stages))
}
result
}bootstrap_hazard_data()
This function randomly samples one sibling per nest from the offspring survival dataset
arguments:
- offspring: the radio-tracked offspring dataset
- num_boot: a unique number for the pbapply simulation (don’t specify)
- species: species name (either “BC” or “WBC”)
- iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
- alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd function.
- max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
- niter: number of iterations specified in the stochastic simulation
# bootstrap_hazard_data()
bootstrap_hazard_data <-
function(offspring, num_boot, species, iter_add, alpha_value, max_time, niter) {
# set attempt to 0 at start of each loop
attempt <- 0
max_time <- max_time
hzd_ss_try_M_lambda <- 1.1
hzd_ss_try_F_lambda <- 1.1
# store simulated estimates only if lambda is less than 1.1
while( (hzd_ss_try_M_lambda > -1.1 | hzd_ss_try_F_lambda > 1.1) && attempt <= 100 ) {
# next attempt
attempt <- attempt + 1
# sample a new offspring dataset containing only one nest member per draw
offspring_boot <-
offspring %>%
group_by(nest_ID) %>%
sample_n(1)
try(
hzd_ss_try_M <- sshzd(Surv(exit, event, entry) ~ exit,
data = filter(offspring_boot, sex == "M"),
alpha = alpha_value),
silent = TRUE
)
try(
hzd_ss_try_F <- sshzd(Surv(exit, event, entry) ~ exit,
data = filter(offspring_boot, sex == "F"),
alpha = alpha_value),
silent = TRUE
)
# store calculated peak (i.e., the apex of the polynomial curve)
try(
hzd_ss_try_M_lambda <- hzd_ss_try_M$lambda
)
# store calculated peak (i.e., the apex of the polynomial curve)
try(
hzd_ss_try_F_lambda <- hzd_ss_try_F$lambda
)
}
# store simulated estimates only if the hzdcurve can be estimated
while( (exists("hzd_curve_try_M") == FALSE | exists("hzd_curve_try_F") == FALSE) && attempt <= max_time) {
time_vector <- seq(0, max_time - attempt, 1)
# next attempt
attempt <- attempt + 1
try(
hzd_ss_try_M <- sshzd(Surv(exit, event, entry) ~ exit,
data = filter(offspring_boot, sex == "M"),
alpha = alpha_value),
silent = TRUE
)
try(
hzd_ss_try_F <- sshzd(Surv(exit, event, entry) ~ exit,
data = filter(offspring_boot, sex == "F"),
alpha = alpha_value),
silent = TRUE
)
# simulate an estimate
try(
hzd_curve_try_M <-
hzdcurve.sshzd(object = hzd_ss_try_M, time = time_vector, se = TRUE),
silent = TRUE
)
try(
hzd_curve_try_F <-
hzdcurve.sshzd(object = hzd_ss_try_F, time = time_vector, se = TRUE),
silent = TRUE
)
}
# make a list including elements that will be used in the next function
out <- list(offspring_boot = offspring_boot,
iter = num_boot + ((iter_add - 1) * niter),
species = species,
time_vector = time_vector,
lambdaM = hzd_ss_try_M_lambda,
lambdaF = hzd_ss_try_F_lambda)
}run_bootstrap_juv_hazd_ad_surv_ASR()
initiates the bootstrap_hazard_data() and runs the juvenile survival analysis to derive sex- and stage-specific survival rates. Then constructs the matrix model with the rates and projects this into the future to acquire the stable stage distribution and derive the ASR.
arguments:
- num_boot: a unique number for the pbapply simulation (don’t specify)
- offspring: the radio-tracked offspring dataset
- k_dist: a two-number vector containing the mean and sd of clutch size (output from the analyses below)
- HSR_dist: a two-number vector containing the mean and sd of hatching sex ratio (output from the analyses below)
- h_dist: a two-number vector containing the mean and sd of the harem size index (output from the analyses below)
- egg_surv_dist: a two-number vector containing the mean and sd of egg survival (output from the analyses below)
- flight_age_distF: a two-number vector containing the mean and sd of female-specific flight age (output from the analyses below)
- flight_age_distM: a two-number vector containing the mean and sd of male-specific flight age (output from the analyses below)
- fledge_age_distF: a two-number vector containing the mean and sd of female-specific fledge age (output from the analyses below)
- fledge_age_distM: a two-number vector containing the mean and sd of male-specific fledge age (output from the analyses below)
- bootstrap_name: a unique string for naming of files in the output directory
- adult_surival_boot_out: the bootstrapped adult encounter history data
- species: species name (either “BC” or “WBC”)
- iter_add: incase of parallel simulations on multiple R Sessions, add a unique number here to prevent overwriting output data on disk
- prefix_number: a unique prefix to prevent overwriting output data on disk (e.g., “boot_one_”)
- alpha_value: Parameter defining cross-validation score for smoothing parameter selection in the gss::sshzd function (default is 1.4).
- max_time: maximum number of days since hatch to model survival (‘70’ in the case of our analysis)
- niter: number of iterations specified in the stochastic simulation
- out.directory: location in the project directory where the output files should be stored (e.g., “tmp”)
run_bootstrap_juv_hazd_ad_surv_ASR <- function(num_boot,
offspring,
k_dist,
HSR_dist,
h_dist,
egg_surv_dist,
flight_age_distF,
flight_age_distM,
fledge_age_distF,
fledge_age_distM,
bootstrap_name,
adult_surival_boot_out,
species,
iter_add,
prefix_number,
alpha_value = 1.4,
max_time = 70,
niter,
output.directory)
{
# run the sampling function and specify the datasets
coucal_boot_list <-
bootstrap_hazard_data(offspring = offspring,
num_boot = num_boot,
species = species,
iter_add = iter_add,
alpha_value = alpha_value,
max_time = max_time,
niter = niter)
# run the survival analysis and ASR deduction on the sampled data
# specify the bootstrapped data samples (from the previous function)
offspring_data <- coucal_boot_list[["offspring_boot"]]
# clean up capture histories
offspring_data <-
offspring_data %>%
ungroup() %>%
as.data.frame()
# fit smoothed spline of hazard function for either sex
M_haz_ss <- sshzd(Surv(exit, event, entry) ~ exit,
data = filter(offspring_data, sex == "M"),
alpha = alpha_value)
F_haz_ss <- sshzd(Surv(exit, event, entry) ~ exit,
data = filter(offspring_data, sex == "F"),
alpha = alpha_value)
haz_ss_lambda <- list(Male_haz_ss_lambda = M_haz_ss$lambda,
Female_haz_ss_lambda = F_haz_ss$lambda)
# extract fitted estimates from the spline function
M_haz_ss_curve <-
hzdcurve.sshzd(object = M_haz_ss, time = coucal_boot_list[["time_vector"]], se = TRUE)
F_haz_ss_curve <-
hzdcurve.sshzd(object = F_haz_ss, time = coucal_boot_list[["time_vector"]], se = TRUE)
haz_ss_curve <-
expand.grid(species = species,
age = coucal_boot_list[["time_vector"]],
sex = c("Male", "Female")) %>%
mutate(fit = c(M_haz_ss_curve$fit, F_haz_ss_curve$fit),
se = c(M_haz_ss_curve$se, F_haz_ss_curve$se)) %>%
mutate(estimate = 1 - fit,
upper = 1 - fit * exp(1.96 * se),
lower = 1 - fit / exp(1.96 * se),
iter = num_boot)
# transform the daily nestling survival (DCS) to apparent fledgling success
# by calculating the product of all DCS estimates:
fledge_ageM <- rnorm(1, fledge_age_distM[1], fledge_age_distM[2])
fledge_ageF <- rnorm(1, fledge_age_distF[1], fledge_age_distF[2])
flight_ageM <- rnorm(1, flight_age_distM[1], flight_age_distM[2])
flight_ageF <- rnorm(1, flight_age_distF[1], flight_age_distF[2])
coucal_nestling_survivalF <-
haz_ss_curve %>%
filter(age <= fledge_ageF & sex == "Female") %>%
group_by(sex) %>%
dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
mutate(stage = "nestling",
rate = "survival")
coucal_nestling_survivalM <-
haz_ss_curve %>%
filter(age <= fledge_ageM & sex == "Male") %>%
group_by(sex) %>%
dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
mutate(stage = "nestling",
rate = "survival")
coucal_nestling_survival <-
bind_rows(coucal_nestling_survivalF, coucal_nestling_survivalM)
coucal_groundling_survivalF <-
haz_ss_curve %>%
filter(age < (flight_ageF + fledge_ageF) & age > fledge_ageF & sex == "Female") %>%
group_by(sex) %>%
dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
mutate(stage = "groundling",
rate = "survival")
coucal_groundling_survivalM <-
haz_ss_curve %>%
filter(age < (flight_ageM + fledge_ageM) & age > fledge_ageM & sex == "Male") %>%
group_by(sex) %>%
dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
mutate(stage = "groundling",
rate = "survival")
coucal_groundling_survival <-
bind_rows(coucal_groundling_survivalF, coucal_groundling_survivalM)
coucal_fledgling_survivalF <-
haz_ss_curve %>%
filter(age >= (flight_ageF + fledge_ageF) & sex == "Female") %>%
group_by(sex) %>%
dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
mutate(stage = "fledgling",
rate = "survival")
coucal_fledgling_survivalM <-
haz_ss_curve %>%
filter(age >= (flight_ageM + fledge_ageM) & sex == "Male") %>%
group_by(sex) %>%
dplyr::summarise(value = prod(estimate), .groups = 'drop') %>%
mutate(stage = "fledgling",
rate = "survival")
coucal_fledgling_survival <-
bind_rows(coucal_fledgling_survivalF, coucal_fledgling_survivalM)
coucal_adult_survival <-
filter(adult_surival_boot_out, iter == sample(1:niter, 1)) %>%
dplyr::select(-species, -iter) %>%
mutate(stage = "adult",
rate = "survival")
egg_survival <- rnorm(1, egg_surv_dist[1], egg_surv_dist[2])
coucal_egg_survival <-
data.frame(sex = NA,
value = egg_survival,
stage = c("egg"),
rate = c("survival"))
coucal_mating_system <-
data.frame(sex = NA,
value = 1/rnorm(1, h_dist[1], h_dist[2]),
stage = c("h"),
rate = c("fecundity"))
coucal_clutch_size <-
data.frame(sex = NA,
value = rnorm(1, k_dist[1], k_dist[2]),
stage = c("k"),
rate = c("fecundity"))
coucal_HSR <-
data.frame(sex = NA,
value = rnorm(1, HSR_dist[1], HSR_dist[2]),
stage = c("HSR"),
rate = c("fecundity"))
coucal_flight_age <-
data.frame(sex = c("Female", "Male"),
value = c(flight_ageF, flight_ageM),
stage = c("flight_age"),
rate = c("development"))
coucal_fledge_age <-
data.frame(sex = c("Female", "Male"),
value = c(fledge_ageF, fledge_ageM),
stage = c("fledge_age"),
rate = c("development"))
# Bind the juvenile and adult dataframe with the nestlings
coucal_vital_rates <-
bind_rows(coucal_egg_survival,
coucal_nestling_survival,
coucal_groundling_survival,
coucal_fledgling_survival,
coucal_adult_survival,
coucal_mating_system,
coucal_clutch_size,
coucal_HSR,
coucal_flight_age,
coucal_fledge_age) %>%
as.data.frame() %>%
mutate(iter = num_boot, #+ ((iter_add - 1) * niter),
species = species) %>%
arrange(sex, rate, stage)
# Create a list of demographic rates from the survival analyses above
demographic_rates <- list(Egg_survival = pull(filter(coucal_vital_rates, stage == "egg"), value),
F_Nestling_survival = pull(filter(coucal_vital_rates,
stage == "nestling" & rate == "survival" & sex == "Female"),
value),
F_Groundling_survival = pull(filter(coucal_vital_rates,
stage == "groundling" & rate == "survival" & sex == "Female"),
value),
F_Fledgling_survival = pull(filter(coucal_vital_rates,
stage == "fledgling" & rate == "survival" & sex == "Female"),
value),
F_Adult_survival = pull(filter(coucal_vital_rates,
stage == "adult" & rate == "survival" & sex == "Female"),
value),
M_Nestling_survival = pull(filter(coucal_vital_rates,
stage == "nestling" & rate == "survival" & sex == "Male"),
value),
M_Groundling_survival = pull(filter(coucal_vital_rates,
stage == "groundling" & rate == "survival" & sex == "Male"),
value),
M_Fledgling_survival = pull(filter(coucal_vital_rates,
stage == "fledgling" & rate == "survival" & sex == "Male"),
value),
M_Adult_survival = pull(filter(coucal_vital_rates,
stage == "adult" & rate == "survival" & sex == "Male"),
value),
# Define hatching sex ratio
HSR = pull(filter(coucal_vital_rates, stage == "HSR"), value),
# Define the mating system (h), and clutch size (k)
h = pull(filter(coucal_vital_rates, stage == "h"), value),
k = pull(filter(coucal_vital_rates, stage == "k"), value))
# Build matrix based on rates specified in the list above
demographic_matrix <- coucal_matrix(demographic_rates_ = demographic_rates)
# Determine the ASR at the stable stage distribution
ASR_SSD <- matrix_ASR(M = demographic_matrix,
h = demographic_rates$h,
HSR = demographic_rates$HSR,
iterations = 100,
num_boot = num_boot,
species = species,
iter_add = 1,
n = rep(10, nrow(demographic_matrix)))
# Extract ASR
ASR_estimate <- ASR_SSD$ASR
# make a list of all the results from this iteration
bootstrap_results_list <-
list(offspring_hazard_lambda = haz_ss_lambda,
offspring_hazard_rates = haz_ss_curve,
coucal_vital_rates = coucal_vital_rates,
ASR_SSD = ASR_SSD,
bootstrapped_data = coucal_boot_list)
# extract file directory
file_directory <- paste0(getwd(), output.directory, "/", bootstrap_name, "_", num_boot, ".Rds") # DEMO "/output/bootstraps/hazard/raw/"
# save the results of the iteration
save(bootstrap_results_list, file = file_directory)
result = bootstrap_results_list
result
}hazard_boot_out_wrangle()
If the product of run_bootstrap_juv_hazd_ad_surv_ASR() is saved on disk, this function will tidy up the components of the raw output for subsequent analyses
arguments:
- species: species name (either “BC” or “WBC”)
- niter: number of iterations specified in the stochastic simulation
- output_dir: location in the project directory where the output files are to be found by the function (e.g., “tmp”)
- rds_file: name of the stochastic simulation rds file saved to disk
# single_model_boot_out_wrangle()
hazard_boot_out_wrangle <-
function(species, niter, output_dir, rds_file){
# specify directory string
boot_path <- paste0(output_dir, species, rds_file, ".rds")
# load each of the four bootstrap output files
bootstrap_out <- readRDS(file = boot_path)
# bind all the fledgling daily survival rates output together
hazard_rates_boot <-
do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
function(x) bootstrap_out["offspring_hazard_rates", x]$offspring_hazard_rates)) %>%
arrange(iter, sex, age)
# bind all the ASR estimates output together
ASR_ests_boot <-
do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
function(x) bootstrap_out["ASR_SSD", x]$ASR_SSD$ASR)) %>%
as.data.frame() %>%
mutate(species = species)
# bind all the ASR estimates output together
vital_rate_ests_boot <-
do.call(rbind, lapply(seq(from = 1, to = niter, by = 1),
function(x) bootstrap_out["coucal_vital_rates", x]$coucal_vital_rates))
# tidy all the bound data together into a mega list of results
boot_out_tidy <-
list(hazard_rates_boot = hazard_rates_boot,
vital_rate_ests_boot = vital_rate_ests_boot,
ASR_ests_boot = ASR_ests_boot)
boot_out_tidy
}sex_diff_hazard()
takes the list of results from the bootstrap procedure and calculates the sex difference in survival at each life stage for each iteration
arguments:
- boot_out_list: the object produced by the hazard_boot_out_wrangle() above
- niter: number of iterations specified in the stochastic simulation
sex_diff_hazard <- function(boot_out_list, niter) {
# make an empty datarame to store the results
sex_diff_surv_output <- data.frame(Adult = niter,
Fledgling = niter,
Groundling = niter,
Nestling = niter,
ISR = niter,
HSR = niter,
h = niter)
# for loop to go through each iteration and calculate the difference between
# female and male survival rates for each stage.
for(i in 1:niter){
Adult <-
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "adult" & rate == "survival" & sex == "Female"), value) -
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "adult" & rate == "survival" & sex == "Male"), value)
Fledgling <-
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "fledgling" & rate == "survival" & sex == "Female"), value) -
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "fledgling" & rate == "survival" & sex == "Male"), value)
Groundling <-
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "groundling" & rate == "survival" & sex == "Female"), value) -
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "groundling" & rate == "survival" & sex == "Male"), value)
Nestling <-
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "nestling" & rate == "survival" & sex == "Female"), value) -
pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "nestling" & rate == "survival" & sex == "Male"), value)
HSR <-
0.5 - pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "HSR"), value)
h <-
1 - pull(filter(boot_out_list$vital_rate_ests_boot[which(boot_out_list$vital_rate_ests_boot$iter == i), ],
stage == "h"), value)
sex_diff_surv_output[i, 1] <- ifelse(length(Adult) > 0, Adult, NA)
sex_diff_surv_output[i, 2] <- ifelse(length(Fledgling) > 0, Fledgling, NA)
sex_diff_surv_output[i, 3] <- ifelse(length(Groundling) > 0, Groundling, NA)
sex_diff_surv_output[i, 4] <- ifelse(length(Nestling) > 0, Nestling, NA)
sex_diff_surv_output[i, 6] <- ifelse(length(HSR) > 0, HSR, NA)
sex_diff_surv_output[i, 7] <- ifelse(length(h) > 0, h, NA)
}
# restructure the output and lable columns
sex_diff_surv_output <- reshape2::melt(data = sex_diff_surv_output)
colnames(sex_diff_surv_output) <- c("stage", "difference")
# return the output
sex_diff_surv_output
}Wrangling survival data
The following chunks show how we wrangled the raw field data to consolidate the data into the appropriate format for the survival analyses
Offspring
status_dat_all <-
# read raw data
read.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv",
header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
# rename ring_ID column
dplyr::rename(ring_ID = Ring_ID) %>%
# make all entries lower case for consistency
mutate(Fledged_status = tolower(Fledged.),
site = tolower(site)) %>%
# select variables of interest
dplyr::select(species, ring_ID, lab_no, sex, year, site, nest_ID, pref_age,
Fledged_status, postf_age, postf_status, ageC, statusC,
lay_date, hatch_order) %>%
# remove all white space from data
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
# specify empty data as NA
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
# exclude all individuals that died in the nest
# filter(Fledged_status == "yes") %>%
# classify columns
mutate(sex = as.factor(sex),
ageC = as.numeric(ageC),
postf_age = as.numeric(postf_age),
postf_status = as.numeric(postf_status),
hatch_order = as.numeric(hatch_order),
pref_age = as.numeric(pref_age)) %>%
# remove rows with missing sex, age, and status info
filter(!is.na(sex) & !is.na(ageC) & !is.na(statusC)) %>%
# make a unique id for each individual
mutate(ind_ID = paste(nest_ID, lab_no, ring_ID, sep = "_"),
# create the age of entry into the data (all at age 15)
entry = 0,
# specify the age of death or censoring
exit = ageC,
# make the event numeric and specify if
# the individual died (1) or was censored (0)
event = as.numeric(statusC)) %>%
# consolidate to variables of interest
dplyr::select(species, ind_ID, nest_ID, year, sex, entry, exit, event, hatch_order)Adult annual survival
detect_dat <-
read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>%
dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>%
mutate(month = str_sub(date_dec, start = 5, end = 6),
day = str_sub(date_dec, start = 7, end = 8)) %>%
mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>%
dplyr::select(-month, -day, -date_dec) %>%
mutate(across(c("sex", "site", "age_status"), tolower)) %>%
mutate(sex = ifelse(sex == "female", "F", ifelse(sex == "male", "M", "XXX")),
age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>%
mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor))
# assess the first and last year of study for each species
detect_dat %>%
dplyr::group_by(species) %>%
dplyr::summarise(first_year = min(year),
last_year = max(year))# A tibble: 3 × 3
species first_year last_year
<fct> <chr> <chr>
1 BC 2001 2018
2 CTC 2016 2019
3 WBC 2005 2019 # make annual capture histories for adults
BC_detect_dat_A <-
detect_dat %>%
filter(age_status == "A" & species == "BC") %>%
mutate(year = as.numeric(year))
BC_detect_dat_A %>%
group_by(ring_ID) %>%
dplyr::summarise(n_years = n_distinct(year)) %>%
arrange(desc(n_years))# A tibble: 280 × 2
ring_ID n_years
<fct> <int>
1 GN41752 2
2 GN52308 2
3 GN55932 2
4 GN76761 2
5 GN76762 2
6 bar_tailed_female 1
7 BCF#34 1
8 BCF#35 1
9 BCF#37 1
10 BCF#39 1
# ℹ 270 more rows# use the BaSTA function "CensusToCaptHist()" function to convert long format
# encounter histories of each individual, to wide format with 1's and 0's for
# each year of encounter
BC_detect_dat_A_ch <-
CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,
d = BC_detect_dat_A$year) %>%
dplyr::rename(ring_ID = ID) %>%
mutate(ID = rownames(.)) %>%
left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>%
distinct()
Black_Coucal_adult_CJS_ch <-
data.frame(ch = apply(BC_detect_dat_A_ch[, 2:19 ] , 1, paste, collapse = "")) %>%
bind_cols(., dplyr::select(BC_detect_dat_A_ch, sex, ring_ID)) %>%
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
distinct()
WBC_detect_dat_A <-
detect_dat %>%
filter(age_status == "A" & species == "WBC") %>%
mutate(year = as.numeric(year))
WBC_detect_dat_A_ch <-
CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,
d = WBC_detect_dat_A$year) %>%
dplyr::rename(ring_ID = ID) %>%
mutate(ID = rownames(.)) %>%
left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>%
distinct()
White_browed_Coucal_adult_CJS_ch <-
data.frame(ch = apply(WBC_detect_dat_A_ch[, 2:16 ] , 1, paste, collapse = "")) %>%
bind_cols(., dplyr::select(WBC_detect_dat_A_ch, sex, ring_ID)) %>%
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
distinct()Bootstrapped mark-recapture analysis of White-browed Coucal adults
Here we resample the annual encounter history of adult White-browed Coucals and conduct a CJS mark-recapture analysis to derive the mean estimate and variance in annual apparent survival.
Set the ‘niter’ to run the bootstrap over several iterations (e.g., 10000). For reproducibility, set the seed with the set.seed() function. Set the ‘output.directory’ so that the subsequent chunks can refer to the appropriate folder within the project directory.
niter = 10000
set.seed(14)
output.directory = "output/bootstraps/"run bootstrap procedure on White-browed Coucals
WBC_adult_survival_CJS_bootstrap_result <-
pbsapply(1:niter, run_bootstrap_adult_survival_CJS,
adult = White_browed_Coucal_adult_CJS_ch,
species = "WBC",
iter_add = 1,
first_year = 2005,
bootstrap_name = "WBC_boot_one",
prefix_number = "boot_one_",
out_directory = "output/bootstraps/")# save the output
saveRDS(WBC_adult_survival_CJS_bootstrap_result,
file = "output/bootstraps/WBC_adult_survival_CJS_bootstrap_result.rds")
# tidy it up
WBC_adult_survival_CJS_bootstrap_tidy <-
surv_boot_out_wrangle(species = "WBC", niter = niter, output_directory = output.directory)
# save the tidy version
saveRDS(WBC_adult_survival_CJS_bootstrap_tidy,
file = "output/bootstraps/WBC_adult_survival_CJS_bootstrap_tidy.rds")Derive the Black Coucal sex-specific survival rates by transforming White-browed Coucal sex-specific survival rates to reflect the ~0.7 ASR that is observed annually
# bootstrapped WBC survival rates
WBC_adult_surival_boot_out <-
WBC_boot_out$adult_reals_survival_rates_boot %>%
dplyr::rename(value = estimate)
# boostrapped BC survival rates (i.e., derived from the WBC rates)
trans_WBC_adult_surival_boot_out <-
WBC_boot_out$adult_reals_survival_rates_boot %>%
dplyr::rename(value = estimate) %>%
mutate(value = ifelse(sex == "Female", value * 0.6, value * 1.4))Parameterizing demographic rates
Each of the following rates has a mean and standard deviation that are used in the stochastic simulation of adult sex ratio.
Mating system
mating_dat <-
# read raw data
read.csv("data/raw/Coucal_Nr_mates_2001_2019_20200123.csv",
header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
# make all entries lower case for consistency
mutate(site = tolower(site),
age_status = tolower(age_status),
sex = tolower(sex)) %>%
# specify empty data as NA
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
# classify columns
mutate(sex = as.factor(sex),
species = as.factor(species),
Nr_partners = as.numeric(Nr_partners),
age_status = as.factor(age_status),
site = as.factor(site)) %>%
# remove rows with missing sex and age, and post-fledged status info
filter(!is.na(sex) & !is.na(Nr_partners) & !is.na(species) & species != "CTC") %>%
# make a unique id for each individual
mutate(ind_ID = str_replace_all(ring_ID, fixed(" "), "")) %>%
mutate(ID_partner.s. = str_replace_all(ID_partner.s., fixed(" "), "")) %>%
# make a unique id for each individual
mutate(ind_ID = str_replace_all(ind_ID, fixed("_"), "")) %>%
mutate(ID_partner.s. = str_replace_all(ID_partner.s., fixed("_"), "")) %>%
# consolidate to variables of interest
dplyr::select(year, species, ind_ID, site, sex, age_status, Nr_partners, ID_partner.s.)
# sex_specific_mating_system <-
# mating_dat %>%
# group_by(species, sex) %>%
# dplyr::summarise(mean_annual_no_mates = mean(Nr_partners, na.rm = TRUE),
# var_annual_no_mates = var(Nr_partners, na.rm = TRUE),
# median_annual_no_mates = median(Nr_partners, na.rm = TRUE),
# sd_annual_no_mates = sd(Nr_partners, na.rm = TRUE),
# n = n_distinct(ind_ID), .groups = "drop",
# n_years = n_distinct(year)) %>%
# mutate(sample_size = paste("n = ", n, sep = ""),
# species_plot = ifelse(species == "WBC", 1.9, 1.1),
# species_lab = ifelse(species == "WBC", 1, 2))
sex_specific_mating_system <-
mating_dat %>%
group_by(ind_ID, sex, species) %>%
dplyr::summarise(mean_annual_no_mates_ind = mean(Nr_partners, na.rm = TRUE),
n_years = n_distinct(year)) %>%
mutate(mu_i = ifelse(mean_annual_no_mates_ind <= 1, 1, mean_annual_no_mates_ind)) %>%
group_by(species, sex) %>%
dplyr::summarise(mean_annual_no_mates = mean(mu_i, na.rm = TRUE),
var_annual_no_mates = var(mu_i, na.rm = TRUE),
median_annual_no_mates = median(mu_i, na.rm = TRUE),
sd_annual_no_mates = sd(mu_i, na.rm = TRUE),
n = n_distinct(ind_ID), .groups = "drop") %>%
mutate(sample_size = paste("n = ", n, sep = ""),
species_plot = ifelse(species == "WBC", 1.9, 1.1),
species_lab = ifelse(species == "WBC", 1, 2))
# To obtain a female-based h index for each population the inverse of the mean
# mu is calculated (i.e., Eq. 2 in manuscript)
BC_h <-
as.numeric(sex_specific_mating_system[which(
sex_specific_mating_system$sex == "female" &
sex_specific_mating_system$species == "BC"), "mean_annual_no_mates"])
WBC_h <-
as.numeric(sex_specific_mating_system[which(
sex_specific_mating_system$sex == "female" &
sex_specific_mating_system$species == "WBC"), "mean_annual_no_mates"])
BC_h_sd <-
as.numeric(sex_specific_mating_system[which(
sex_specific_mating_system$sex == "female" &
sex_specific_mating_system$species == "BC"), "sd_annual_no_mates"])
WBC_h_sd <-
as.numeric(sex_specific_mating_system[which(
sex_specific_mating_system$sex == "female" &
sex_specific_mating_system$species == "WBC"), "sd_annual_no_mates"])
coucal_mating_system <-
data.frame(trait = c("mating_system"),
species = c("BC", "WBC"),
mean = c(BC_h, WBC_h),
sd = c(BC_h_sd, WBC_h_sd))
mating_dat$species <-
factor(mating_dat$species ,
levels = c("BC",
"WBC"))
mating_dat_plotting <-
mating_dat %>%
mutate(species_dummy = ifelse(species == "WBC", 1, 0)) %>%
mutate(species_plot = ifelse(species == "WBC", 2.1, 0.9))
mating_system_plot <-
ggplot2::ggplot() +
theme_bw() +
geom_jitter(aes(y = Nr_partners, x = species_plot,
color = sex),
data = mating_dat_plotting,
width = 0.1, height = 0.1,
alpha = 0.75,
fill = "white", shape = 16) +
geom_pointrange(data = sex_specific_mating_system,
aes(y = mean_annual_no_mates, x = species_plot,
ymin = (mean_annual_no_mates-sd_annual_no_mates),
ymax = (mean_annual_no_mates+sd_annual_no_mates),
fill = sex),
size = 0.8, shape = 21, color = "black") +
geom_text(aes(y = 0.2, x = species_lab, label = sample_size),
data = sex_specific_mating_system,
size = 3, family = "Verdana") +
geom_hline(yintercept = 1.5, linetype = "dashed",
alpha = 0.5, color = "grey20") +
facet_grid(. ~ sex, labeller = as_labeller(sex_names)) +
luke_theme +
theme(legend.position = "none",
axis.title.x = element_blank(),
strip.background = element_blank(),
strip.text = element_text(size = 12, face = "italic")) +
scale_x_continuous(limits = c(0.5, 2.5),
breaks = c(1, 2), labels = species_names) +
scale_y_continuous(limits = c(0, 5.5), expand = c(0, 0),
breaks = c(1, 2, 3, 4, 5),
labels = c(1, 2, 3, 4, 5)) +
ylab("Number of unique mates\nper year (± 1 SD)") +
scale_color_manual(values = plot_palette_sex) +
scale_fill_manual(values = plot_palette_sex) +
annotate(geom = "text", y = 1.25, x = 1.5,
label = "monogamy", family = "Verdana",
color = "black", size = 3, fontface = 'italic', hjust = 0.5) +
annotate(geom = "text", y = 1.75, x = 1.5,
label = "polygamy", family = "Verdana",
color = "black", size = 3, fontface = 'italic', hjust = 0.5)
Clutch size
# import raw csv data into R
egg_data <-
read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls", na = "NA", col_types = "text") %>%
dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks,
nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
`unkn_sex_eggs+chiks`) %>%
filter(species != "CTC") %>%
filter(!is.na(nest_success)) %>%
dplyr::rename(n_eggs = `Total clutch size`,
n_chicks = all_chicks,
male_eggs = `male_chi+eggs`,
female_eggs = `female_chi+eggs`,
unknown_eggs = `unkn_sex_eggs+chiks`) %>%
dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
n_chicks = na_if(n_chicks, "?")) %>%
dplyr::mutate(n_eggs = as.numeric(n_eggs),
n_chicks = as.numeric(n_chicks),
male_eggs = as.numeric(male_eggs),
female_eggs = as.numeric(female_eggs),
unknown_eggs = as.numeric(unknown_eggs)) %>%
filter(n_eggs >= n_chicks)
#### clutch size summary ----
coucal_clutch_size <-
data.frame(trait = c("clutch_size"),
species = c("BC", "WBC"),
mean = c(filter(egg_data, species == "BC") %>%
summarise(mean_n_eggs = mean(n_eggs, na.rm = TRUE)) %>%
pull(mean_n_eggs),
filter(egg_data, species == "WBC") %>%
summarise(mean_n_eggs = mean(n_eggs, na.rm = TRUE)) %>%
pull(mean_n_eggs)),
sd = c(filter(egg_data, species == "BC") %>%
summarise(mean_n_eggs = sd(n_eggs, na.rm = TRUE)) %>%
pull(mean_n_eggs),
filter(egg_data, species == "WBC") %>%
summarise(mean_n_eggs = sd(n_eggs, na.rm = TRUE)) %>%
pull(mean_n_eggs)),
n_nests = c(filter(egg_data, species == "BC") %>%
summarise(n_nests = n_distinct(nest_ID)) %>%
pull(n_nests),
filter(egg_data, species == "WBC") %>%
summarise(n_nests = n_distinct(nest_ID)) %>%
pull(n_nests)),
n_years = c(filter(egg_data, species == "BC") %>%
summarise(n_nests = n_distinct(year)) %>%
pull(n_nests),
filter(egg_data, species == "WBC") %>%
summarise(n_nests = n_distinct(year)) %>%
pull(n_nests)))Egg hatchability
# import raw csv data into R
egg_data <-
read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls",
na = "NA", col_types = "text") %>%
dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks,
nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
`unkn_sex_eggs+chiks`) %>%
filter(species != "CTC") %>%
filter(!is.na(nest_success)) %>%
dplyr::rename(n_eggs = `Total clutch size`,
n_chicks = all_chicks,
male_eggs = `male_chi+eggs`,
female_eggs = `female_chi+eggs`,
unknown_eggs = `unkn_sex_eggs+chiks`) %>%
dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
n_chicks = na_if(n_chicks, "?")) %>%
dplyr::mutate(n_eggs = as.numeric(n_eggs),
n_chicks = as.numeric(n_chicks),
male_eggs = as.numeric(male_eggs),
female_eggs = as.numeric(female_eggs),
unknown_eggs = as.numeric(unknown_eggs)) %>%
filter(n_eggs >= n_chicks) %>%
arrange(nest_ID)
#### egg hatchability summary ----
egg_surv_data <-
egg_data %>%
dplyr::mutate(n_dead_eggs = n_eggs - n_chicks,
n_alive_eggs = n_chicks) %>%
arrange(nest_ID)
#### sample size summary ----
egg_surv_data %>%
dplyr::group_by(species) %>%
dplyr::summarise(n_nests = n_distinct(nest_ID))# A tibble: 2 × 2
species n_nests
<chr> <int>
1 BC 271
2 WBC 213BC_egg_survival <-
lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~
1 +
(1| nest_ID),
data = filter(egg_surv_data, species == "BC"),
family = binomial)
WBC_egg_survival <-
lme4::glmer(cbind(n_alive_eggs, n_dead_eggs) ~
1 +
(1| nest_ID),
data = filter(egg_surv_data, species == "WBC"),
family = binomial)
# # run tidy bootstrap to obtain model diagnostics
tidy_BC_egg_survival <-
tidy(BC_egg_survival, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000
tidy_WBC_egg_survival <-
tidy(WBC_egg_survival, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000
coucal_egg_survival <-
data.frame(trait = c("egg_survival"),
species = c("BC", "WBC"),
mean = c(invlogit(model_parameters(BC_egg_survival)$Coefficient)[1],
invlogit(model_parameters(WBC_egg_survival)$Coefficient)[1]),
CI_low = c(invlogit(model_parameters(BC_egg_survival)$CI_low)[1],
invlogit(model_parameters(WBC_egg_survival)$CI_low)[1]),
CI_high = c(invlogit(model_parameters(BC_egg_survival)$CI_high)[1],
invlogit(model_parameters(WBC_egg_survival)$CI_high)[1]),
n_nests = c(filter(egg_data, species == "BC") %>%
summarise(n_nests = n_distinct(nest_ID)) %>%
pull(n_nests),
filter(egg_data, species == "WBC") %>%
summarise(n_nests = n_distinct(nest_ID)) %>%
pull(n_nests)),
n_years = c(filter(egg_data, species == "BC") %>%
summarise(n_nests = n_distinct(year)) %>%
pull(n_nests),
filter(egg_data, species == "WBC") %>%
summarise(n_nests = n_distinct(year)) %>%
pull(n_nests))) %>%
mutate(sd = ifelse(!is.na(CI_low),
approx_sd(x1 = CI_low, x2 = CI_high),
CI_low))Hatching sex ratio
egg_data <-
read_xls("data/raw/Egg_surv_data_2001_2019_20210524.xls",
na = "NA", col_types = "text") %>%
dplyr::select(species, nest_ID, year, `Total clutch size`, all_chicks,
nst_suc_sho, nest_success, `male_chi+eggs`, `female_chi+eggs`,
`unkn_sex_eggs+chiks`) %>%
filter(species != "CTC") %>%
filter(!is.na(nest_success)) %>%
dplyr::rename(n_eggs = `Total clutch size`,
n_chicks = all_chicks,
male_eggs = `male_chi+eggs`,
female_eggs = `female_chi+eggs`,
unknown_eggs = `unkn_sex_eggs+chiks`) %>%
dplyr::mutate(n_eggs = na_if(n_eggs, "?"),
n_chicks = na_if(n_chicks, "?")) %>%
dplyr::mutate(n_eggs = as.numeric(n_eggs),
n_chicks = as.numeric(n_chicks),
male_eggs = as.numeric(male_eggs),
female_eggs = as.numeric(female_eggs),
unknown_eggs = as.numeric(unknown_eggs)) %>%
filter(n_eggs >= n_chicks)
#### hatching sex ratio summary ----
egg_HSR_data <-
egg_data %>%
dplyr::filter(n_eggs == (male_eggs + female_eggs + unknown_eggs)) %>%
dplyr::filter(unknown_eggs == 0) %>%
mutate(HSR = male_eggs/(male_eggs + female_eggs + unknown_eggs)) %>%
mutate(species_plot = ifelse(species == "WBC", 2.2, 0.8),
species_plot2 = ifelse(species == "WBC", 2, 1))
egg_HSR_data %>%
dplyr::group_by(species) %>%
dplyr::summarise(n_distinct(nest_ID))# A tibble: 2 × 2
species `n_distinct(nest_ID)`
<chr> <int>
1 BC 70
2 WBC 51mod_HSR_BC <-
lme4::glmer(cbind(male_eggs, female_eggs) ~
1 +
(1| nest_ID),
data = filter(egg_HSR_data, species == "BC"),
family = binomial)
mod_HSR_WBC <-
lme4::glmer(cbind(male_eggs, female_eggs) ~
1 +
(1| nest_ID),
data = filter(egg_HSR_data, species == "WBC"),
family = binomial)
# # run tidy bootstrap to obtain model diagnostics
tidy_mod_HSR_BC <-
tidy(mod_HSR_BC, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000
tidy_mod_HSR_WBC <-
tidy(mod_HSR_WBC, conf.int = TRUE, conf.method = "boot", nsim = 1000) # DEMO 1000
coucal_HSR <-
data.frame(trait = "hatching_sex_ratio",
species = c("BC", "WBC"),
mean = c(invlogit(model_parameters(mod_HSR_BC)$Coefficient)[1],
invlogit(model_parameters(mod_HSR_WBC)$Coefficient)[1]),
CI_low = c(invlogit(model_parameters(mod_HSR_BC)$CI_low)[1],
invlogit(model_parameters(mod_HSR_WBC)$CI_low)[1]),
CI_high = c(invlogit(model_parameters(mod_HSR_BC)$CI_high)[1],
invlogit(model_parameters(mod_HSR_WBC)$CI_high)[1]),
n_nests = c(filter(egg_HSR_data, species == "BC") %>%
summarise(n_nests = n_distinct(nest_ID)) %>%
pull(n_nests),
filter(egg_HSR_data, species == "WBC") %>%
summarise(n_nests = n_distinct(nest_ID)) %>%
pull(n_nests)),
n_years = c(filter(egg_HSR_data, species == "BC") %>%
summarise(n_nests = n_distinct(year)) %>%
pull(n_nests),
filter(egg_HSR_data, species == "WBC") %>%
summarise(n_nests = n_distinct(year)) %>%
pull(n_nests))) %>%
mutate(sd = ifelse(!is.na(CI_low),
approx_sd(x1 = CI_low, x2 = CI_high),
CI_low),
species_plot = ifelse(species == "WBC", 1.8, 1.2))
HSR_plot <-
ggplot2::ggplot() +
geom_boxplot(data = egg_HSR_data,
aes(x = species_plot, y = HSR,
group = species, fill = species),
color = "grey50",
width = 0.05, alpha = 0.5,
position = position_dodge(width = 0)) +
geom_errorbar(data = coucal_HSR,
aes(x = species_plot, ymax = CI_high, ymin = CI_low),
alpha = 1, color = "black", width = 0.05, lwd = 0.5) +
geom_point(data = coucal_HSR,
aes(x = species_plot, y = mean, fill = species),
lwd = 3, shape = 21, color= "black") +
geom_jitter(data = egg_HSR_data,
aes(x = species_plot2, y = HSR,
group = species,
fill = species, color = species),
width = 0.02, alpha = 0.2, shape = 19) +
luke_theme +
theme(legend.position = "none",
panel.border = element_blank(),
strip.background = element_blank(),
axis.title.x = element_blank(),
axis.text.x = element_text(12),
legend.background = element_blank(),
panel.grid.major.y = element_line(colour = "grey70", size = 0.1),
plot.title = element_text(hjust = 0.5, face = "italic")) +
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.25)) +
scale_x_discrete(labels = c("1" = "Black Coucal",
"2" = "White-browed Coucal")) +
ylab(expression(paste("Hatching sex ratio (prop. male)" %+-% "95% CI", sep = ""))) +
xlab("Species") +
scale_color_manual(values = c(brewer.pal(6, "Dark2")[1], brewer.pal(6, "Set1")[1])) +
scale_fill_manual(values = c(brewer.pal(6, "Dark2")[1], brewer.pal(6, "Set1")[1])) +
ggtitle("Black Coucal White-browed\n Coucal")
Age at fledging (transition from nestling to groundling)
fledge_dat <-
read.csv("data/raw/Coucal_chick_survival_2001-2019_20200129.csv",
header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
dplyr::select(species, sex, Fledge_age, nest_ID, year, lab_no, Ring_ID, site) %>%
filter(site == "Kapunga" & species != "CTC") %>%
mutate(Nst_No = str_replace_all(string = nest_ID, fixed(" "), ""),
Ring_ID = str_replace_all(string = Ring_ID, fixed(" "), ""),
lab_no = str_replace_all(string = lab_no, fixed(" "), "")) %>%
filter(!is.na(Fledge_age)) %>%
mutate(Ind_ID = paste(Nst_No, Ring_ID, lab_no, sep = "_")) %>%
mutate(sex_plot = ifelse(sex == "M", 2.2, 0.8))
#### Modeling (BC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_BC <-
lmer(Fledge_age ~ sex +
(1 | nest_ID) + (1 | year),
data = filter(fledge_dat, species == "BC"))
new_data_BC <-
expand.grid(sex = c("F","M"))
predict_BC <-
predict(mod_fledge_age_BC,
newdata = new_data_BC,
re.form = NA,
se.fit = TRUE,
nsim = 1000)
#### Modeling (WBC) ----
# "Fledge_age" as dependent variable, interaction with sex and species
mod_fledge_age_WBC <-
lmer(Fledge_age ~ sex +
(1 | nest_ID) + (1 | year),
data = filter(fledge_dat, species == "WBC"))
new_data_WBC <-
expand.grid(sex = c("F","M"))
predict_WBC <-
predict(mod_fledge_age_WBC,
newdata = new_data_WBC,
re.form = NA,
se.fit = TRUE,
nsim = 1000)
# tidy up estimates
coucal_fledge_age <-
data.frame(trait = c("fledge_age"),
species = c("BC", "BC", "WBC", "WBC"),
sex = c("F", "M", "F", "M"),
mean = c(predict_BC$fit[1],
predict_BC$fit[2],
predict_WBC$fit[1],
predict_WBC$fit[2]),
CI_low = c(predict_BC$ci.fit[1, 1],
predict_BC$ci.fit[1, 2],
predict_WBC$ci.fit[1, 1],
predict_WBC$ci.fit[1, 2]),
CI_high = c(predict_BC$ci.fit[2, 1],
predict_BC$ci.fit[2, 2],
predict_WBC$ci.fit[2, 1],
predict_WBC$ci.fit[2, 2]),
n_inds = c(filter(fledge_dat, species == "BC" & sex == "F") %>%
summarise(n_ = n_distinct(Ind_ID)) %>%
pull(n_),
filter(fledge_dat, species == "BC" & sex == "M") %>%
summarise(n_ = n_distinct(Ind_ID)) %>%
pull(n_),
filter(fledge_dat, species == "WBC" & sex == "F") %>%
summarise(n_ = n_distinct(Ind_ID)) %>%
pull(n_),
filter(fledge_dat, species == "BC" & sex == "M") %>%
summarise(n_ = n_distinct(Ind_ID)) %>%
pull(n_)),
n_nests = c(filter(fledge_dat, species == "BC" & sex == "F") %>%
summarise(n_ = n_distinct(Nst_No)) %>%
pull(n_),
filter(fledge_dat, species == "BC" & sex == "M") %>%
summarise(n_ = n_distinct(Nst_No)) %>%
pull(n_),
filter(fledge_dat, species == "WBC" & sex == "F") %>%
summarise(n_ = n_distinct(Nst_No)) %>%
pull(n_),
filter(fledge_dat, species == "WBC" & sex == "M") %>%
summarise(n_ = n_distinct(Nst_No)) %>%
pull(n_)),
n_years = c(filter(fledge_dat, species == "BC" & sex == "F") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_),
filter(fledge_dat, species == "BC" & sex == "M") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_),
filter(fledge_dat, species == "WBC" & sex == "F") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_),
filter(fledge_dat, species == "WBC" & sex == "M") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_))) %>%
mutate(sd = ifelse(!is.na(CI_low),
approx_sd(x1 = CI_low, x2 = CI_high),
CI_low))Age at first flight (transition from groundling to fledgling)
flight_dat <-
read.csv("data/raw/coucals_first_flights_data_per_individual2020.csv",
header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
mutate(Nst_No = str_replace_all(string = Nst_No, fixed(" "), ""),
Ind_ID = str_replace_all(string = Ind_ID, fixed(" "), "")) %>%
dplyr::rename(sex = Sex,
nest_ID = Nst_No,
year = Year,
species = Spp,
flight_age = days_since_fledging,
Ring_ID = Ind_ID) %>%
mutate(sex_plot = ifelse(sex == "Male", 2.2, 0.8))
#### Modeling (WBC) ----
# "flight_age" as dependent variable,
# "sex", "fledge_age", and "fledge_mass" as independent
mod_flight_age_WBC <-
lmer(flight_age ~ sex + Fledge_age + Fledge_mass +
(1 | nest_ID) + (1 | year),
data = filter(flight_dat, species == "WBC"))
new_data_WBC <-
expand.grid(sex = c("Female","Male"),
Fledge_age = mean(flight_dat[flight_dat$species == "WBC",]$Fledge_age),
Fledge_mass = mean(flight_dat[flight_dat$species == "WBC",]$Fledge_mass))
predict_WBC <-
predict(mod_flight_age_WBC,
newdata = new_data_WBC,
re.form = NA,
se.fit = TRUE,
nsim = 1000)
#### Modeling (BC) ----
# "flight_age" as dependent variable,
# "sex", "fledge_age", and "fledge_mass" as independent
mod_flight_age_BC <-
lmer(flight_age ~ sex + Fledge_age + Fledge_mass +
(1 | nest_ID) + (1 | year),
data = filter(flight_dat, species == "BC"))
new_data_BC <-
expand.grid(sex = c("Female","Male"),
Fledge_age = mean(flight_dat[flight_dat$species == "BC",]$Fledge_age),
Fledge_mass = mean(flight_dat[flight_dat$species == "BC",]$Fledge_mass))
predict_BC <-
predict(mod_flight_age_BC,
newdata = new_data_BC,
re.form = NA,
se.fit = TRUE,
nsim = 1000)
coucal_flight_age <-
data.frame(trait = c("flight_age"),
species = c("BC", "BC", "WBC", "WBC"),
sex = c("F", "M", "F", "M"),
mean = c(predict_BC$fit[1],
predict_BC$fit[2],
predict_WBC$fit[1],
predict_WBC$fit[2]),
CI_low = c(predict_BC$ci.fit[1, 1],
predict_BC$ci.fit[1, 2],
predict_WBC$ci.fit[1, 1],
predict_WBC$ci.fit[1, 2]),
CI_high = c(predict_BC$ci.fit[2, 1],
predict_BC$ci.fit[2, 2],
predict_WBC$ci.fit[2, 1],
predict_WBC$ci.fit[2, 2]),
n_inds = c(filter(flight_dat, species == "BC" & sex == "Female") %>%
summarise(n_ = n_distinct(Ring_ID)) %>%
pull(n_),
filter(flight_dat, species == "BC" & sex == "Male") %>%
summarise(n_ = n_distinct(Ring_ID)) %>%
pull(n_),
filter(flight_dat, species == "WBC" & sex == "Female") %>%
summarise(n_ = n_distinct(Ring_ID)) %>%
pull(n_),
filter(flight_dat, species == "BC" & sex == "Male") %>%
summarise(n_ = n_distinct(Ring_ID)) %>%
pull(n_)),
n_nests = c(filter(flight_dat, species == "BC" & sex == "Female") %>%
summarise(n_ = n_distinct(nest_ID)) %>%
pull(n_),
filter(flight_dat, species == "BC" & sex == "Male") %>%
summarise(n_ = n_distinct(nest_ID)) %>%
pull(n_),
filter(flight_dat, species == "WBC" & sex == "Female") %>%
summarise(n_ = n_distinct(nest_ID)) %>%
pull(n_),
filter(flight_dat, species == "WBC" & sex == "Male") %>%
summarise(n_ = n_distinct(nest_ID)) %>%
pull(n_)),
n_years = c(filter(flight_dat, species == "BC" & sex == "Female") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_),
filter(flight_dat, species == "BC" & sex == "Male") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_),
filter(flight_dat, species == "WBC" & sex == "Female") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_),
filter(flight_dat, species == "WBC" & sex == "Male") %>%
summarise(n_ = n_distinct(year)) %>%
pull(n_))) %>%
mutate(sd = ifelse(!is.na(CI_low),
approx_sd(x1 = CI_low, x2 = CI_high),
CI_low))Bind all parameters together
parameter_distributions <-
bind_rows(coucal_egg_survival,
coucal_clutch_size,
coucal_HSR,
coucal_mating_system,
coucal_fledge_age,
coucal_flight_age) %>%
dplyr::select(trait, species, sex, mean, sd)Stochastic Demographic Simulation
niter = 10000
set.seed(14)Black Coucals
pull species- and sex-specific means and sds for each parameter
k_dist_BC = c(pull(filter(coucal_clutch_size, species == "BC"), mean),
pull(filter(coucal_clutch_size, species == "BC"), sd))
HSR_dist_BC = c(pull(filter(coucal_HSR, species == "BC"), mean),
pull(filter(coucal_HSR, species == "BC"), sd))
h_dist_BC = c(pull(filter(coucal_mating_system, species == "BC"), mean),
pull(filter(coucal_mating_system, species == "BC"), sd))
egg_surv_dist_BC = c(pull(filter(coucal_egg_survival, species == "BC"), mean),
pull(filter(coucal_egg_survival, species == "BC"), sd))
fledge_age_distF_BC = c(pull(filter(coucal_fledge_age, species == "BC" & sex == "F"),
mean),
pull(filter(coucal_fledge_age, species == "BC" & sex == "F"),
sd))
fledge_age_distM_BC = c(pull(filter(coucal_fledge_age, species == "BC" & sex == "M"),
mean),
pull(filter(coucal_fledge_age, species == "BC" & sex == "M"),
sd))
flight_age_distF_BC = c(pull(filter(coucal_flight_age, species == "BC" & sex == "F"),
mean),
pull(filter(coucal_flight_age, species == "BC" & sex == "F"),
sd))
flight_age_distM_BC = c(pull(filter(coucal_flight_age, species == "BC" & sex == "M"),
mean),
pull(filter(coucal_flight_age, species == "BC" & sex == "M"),
sd))run stochastic demographic simulation
BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc <-
pbsapply(1:niter, run_bootstrap_juv_hazd_ad_surv_ASR,
offspring = status_dat_all %>% filter(species == "BC"),
adult_surival_boot_out = trans_WBC_adult_surival_boot_out,
bootstrap_name = "BC_boot_w_trans_WBC_ad_surv_stoc",
species = "BC",
iter_add = 1,
prefix_number = "boot_w_trans_WBC_ad_surv_stoc_no",
max_time = 70,
niter = niter,
alpha_value = 1.4,
k_dist = k_dist_BC,
HSR_dist = HSR_dist_BC,
h_dist = h_dist_BC,
egg_surv_dist = egg_surv_dist_BC,
fledge_age_distF = fledge_age_distF_BC,
fledge_age_distM = fledge_age_distM_BC,
flight_age_distF = flight_age_distF_BC,
flight_age_distM = flight_age_distM_BC,
output.directory = "/output/bootstraps/hazard/raw")White-browed Coucals
pull species- and sex-specific means and sds for each parameter
k_dist_WBC = c(pull(filter(coucal_clutch_size, species == "WBC"), mean),
pull(filter(coucal_clutch_size, species == "WBC"), sd))
HSR_dist_WBC = c(pull(filter(coucal_HSR, species == "WBC"), mean),
pull(filter(coucal_HSR, species == "WBC"), sd))
h_dist_WBC = c(pull(filter(coucal_mating_system, species == "WBC"), mean),
pull(filter(coucal_mating_system, species == "WBC"), sd))
egg_surv_dist_WBC = c(pull(filter(coucal_egg_survival, species == "WBC"), mean),
pull(filter(coucal_egg_survival, species == "WBC"), sd))
fledge_age_distF_WBC = c(pull(filter(coucal_fledge_age, species == "WBC" & sex == "F"),
mean),
pull(filter(coucal_fledge_age, species == "WBC" & sex == "F"),
sd))
fledge_age_distM_WBC = c(pull(filter(coucal_fledge_age, species == "WBC" & sex == "M"),
mean),
pull(filter(coucal_fledge_age, species == "WBC" & sex == "M"),
sd))
flight_age_distF_WBC = c(pull(filter(coucal_flight_age, species == "WBC" & sex == "F"),
mean),
pull(filter(coucal_flight_age, species == "WBC" & sex == "F"),
sd))
flight_age_distM_WBC = c(pull(filter(coucal_flight_age, species == "WBC" & sex == "M"),
mean),
pull(filter(coucal_flight_age, species == "WBC" & sex == "M"),
sd))run stochastic demographic simulation
WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc <-
pbsapply(1:niter, run_bootstrap_juv_hazd_ad_surv_ASR,
offspring = status_dat_all %>% filter(species == "WBC"),
adult_surival_boot_out = WBC_adult_surival_boot_out,
bootstrap_name = "WBC_boot_w_WBC_ad_surv_stoc",
species = "WBC",
iter_add = 1,
prefix_number = "boot_w_WBC_ad_surv_stoc",
max_time = 70,
niter = niter,
alpha_value = 1.4,
k_dist = k_dist_WBC,
HSR_dist = HSR_dist_WBC,
h_dist = h_dist_WBC,
egg_surv_dist = egg_surv_dist_WBC,
fledge_age_distF = fledge_age_distF_WBC,
fledge_age_distM = fledge_age_distM_WBC,
flight_age_distF = flight_age_distF_WBC,
flight_age_distM = flight_age_distM_WBC,
output.directory = "/output/bootstraps/hazard/raw")Save simulation output
saveRDS(object = BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc,
file = "output/bootstraps/hazard/cooked/BC_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc.rds")
saveRDS(object = WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc,
file = "output/bootstraps/hazard/cooked/WBC_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc.rds")clean up the output from the bootstrap procedure and save as rds
BC_hazard_rate_boot_tidy_trans_no_imm <-
hazard_boot_out_wrangle(species = "BC", niter = 1000,
output_dir = "output/bootstraps/hazard/cooked/",
rds_file =
"_hazard_ASR_bootstrap_result_w_trans_WBC_ad_surv_stoc_no_imm")
WBC_hazard_rate_boot_tidy_no_imm <-
hazard_boot_out_wrangle(species = "WBC", niter = 1000,
output_dir = "output/bootstraps/hazard/cooked/",
rds_file =
"_hazard_ASR_bootstrap_result_w_WBC_ad_surv_stoc_no_imm")save model output
saveRDS(object = BC_hazard_rate_boot_tidy_trans_no_imm,
file =
"output/bootstraps/hazard/cooked/BC_haz_sur_ASR_boot_tidy_stoc_trans_no_imm.rds")
saveRDS(object = WBC_hazard_rate_boot_tidy_no_imm,
file =
"output/bootstraps/hazard/cooked/WBC_haz_sur_ASR_boot_tidy_stoc_no_imm.rds")Sex-biased Demographic stages
# calculate the sex differences in stage specific rates
BC_sex_diff_hazard_output <-
sex_diff_hazard(boot_out_list = BC_hazard_rate_boot_tidy,
niter = 1000) %>%
mutate(species = "BC")
WBC_sex_diff_hazard_output <-
sex_diff_hazard(boot_out_list = WBC_hazard_rate_boot_tidy,
niter = 1000) %>%
mutate(species = "WBC")
# consolidate results
sex_diff_survival_output <-
dplyr::bind_rows(BC_sex_diff_hazard_output,
WBC_sex_diff_hazard_output) %>%
dplyr::filter(stage %in% c("HSR", "Nestling", "Groundling",
"Fledgling", "Adult")) %>%
mutate(stage = factor(stage,
levels = c("HSR", "Nestling", "Groundling",
"Fledgling", "Adult")),
species = factor(species,
levels = c("BC", "WBC")),
difference = -difference)
# check that there are no nosensical values
sex_diff_survival_output %>%
filter(difference < -1 & difference > 1)
ggplot(data = filter(sex_diff_survival_output, stage != "h")) +
geom_histogram(aes(difference), binwidth = 0.01) +
facet_grid(stage ~ species) +
scale_x_continuous(limits = c(-0.5, 0.5)) +
geom_vline(xintercept = 0)
# calculate some summary statistics
sex_diff_survival_summary <-
sex_diff_survival_output %>%
dplyr::group_by(stage, species) %>%
dplyr::summarise(avg = mean(difference, na.rm = TRUE),
median = median(difference, na.rm = TRUE),
var = var(difference, na.rm = TRUE),
max = max(difference, na.rm = TRUE),
min = min(difference, na.rm = TRUE),
sd = sd(difference, na.rm = TRUE))
# specify custom color palette to distingush first-year stages
# (i.e. chicks and juveniles) from adults
cbPalette <- c("#D9D9D9", "#D9D9D9", "#D9D9D9", "#D9D9D9", "#A6A6A6")
species_names <-
c('BC' = "Black Coucal",
'WBC' = "White-browed Coucal")
junk0 <-
data.frame(species = c("BC", "WBC"),
stage = "HSR",
difference = c(NA,NA))
sex_diff_survival_output2 <-
bind_rows(sex_diff_survival_output, junk0) %>%
mutate(difference = as.numeric(difference),
stage = factor(stage,
levels = c("HSR", "Nestling", "Groundling",
"Fledgling", "Adult")),
species = factor(species,
levels = c("BC", "WBC")))
# Figure 2a: plot the sex-biases in survival across the three stages
vital_rate_theme <-
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.position = "none",
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.length = unit(0.1, "cm"),
panel.border = element_blank(),
panel.spacing = unit(0.3, "lines"),
strip.background = element_blank()
)
surv_diff_plot <-
ggplot(aes(y = difference, x = stage, fill = stage),
data = slice_sample(group_by(sex_diff_survival_output, species, stage), n = 1000)) +
geom_violin(draw_quantiles = c(0.025, 0.5, 0.975), color = "grey10",
scale = "width", trim = TRUE, adjust = 1, size = 0.25) +
vital_rate_theme +
facet_grid(. ~ species, labeller = as_labeller(species_names)) +
theme(axis.title.x = element_text(size = 7, colour = "black"),
axis.text.x = element_text(size = 6, angle = 45, hjust = 1, vjust = 1,
colour = "black"),
axis.ticks.x = element_line(size = 0.2),
axis.title.y = element_text(size = 7, hjust = 0.5, vjust = -3),
axis.text.y = element_text(size = 6, colour = "black"),
axis.ticks.y = element_line(size = 0.2),
plot.margin = unit(c(0.2, 1.35, 0.42, 0.3), "cm"),
strip.text = element_text(size = 6, colour = "black", face = "italic")
) +
scale_fill_manual(values = cbPalette) +
scale_y_continuous(limits = c(-0.65, 0.65),
breaks = c(-0.5, -0.25, 0, 0.25, 0.5),
expand = c(0, 0), position = "left",
labels = c("-0.50","-0.25",
expression(phantom("-")*"0.00"),
expression(phantom("-")*"0.25"),
expression(phantom("-")*"0.50"))) +
xlab("Life stage") +
ylab("Sex bias\n") +
scale_x_discrete(labels = c(
"HSR" = "hatching\nsex ratio",
"Nestling" = "nestling\nsurvival",
"Groundling" = "groundling\nsurvival",
"Fledgling" = "fledgling\nsurvival",
"Adult" = "adult\nsurvival"))
background <-
ggplot(aes(y = difference, x = stage, fill = stage),
data = sex_diff_survival_output) +
annotate("rect", xmin = -Inf, xmax = Inf,
ymin = -Inf, ymax = 0, alpha = 0.7,
fill = pull(ggthemes_data$wsj$palettes$colors6[3,2])) +
annotate("rect", xmin = -Inf, xmax = Inf,
ymin = 0, ymax = Inf, alpha = 0.7,
fill = pull(ggthemes_data$wsj$palettes$colors6[2,2])) +
annotate("text", x = c(5), y = c(-0.5),
label = c("\u2640"), size = 4,
family = "Menlo",
vjust = c(0), hjust = c(0.5), colour = "grey10") +
annotate("text", x = c(1), y = c(0.5),
label = c("\u2642"), size = 4,
family = "Menlo",
vjust = c(1), hjust = c(0.5), colour = "grey10") +
facet_grid(. ~ species, labeller = as_labeller(species_names)) +
vital_rate_theme +
theme(axis.title.x = element_text(size = 7, colour = "white"),
axis.text.x = element_text(size = 6, angle = 45, hjust = 1,
vjust = 1, colour = "white"),
axis.ticks.x = element_line(size = 0.2, colour = "white"),
axis.title.y = element_text(size = 7, hjust = 0.5, vjust = -1,
colour = "white"),
axis.text.y = element_text(size = 6, colour = "white"),
axis.ticks.y = element_line(size = 0.2, colour = "white"),
plot.margin = unit(c(0.2, 1.35, 1.55, 0.3), "cm"),
strip.text = element_text(size = 6, colour = "white")
) +
scale_x_discrete(labels = c(
"HSR" = "hatching\nsex ratio",
"Nestling" = "nestling\nsurvival",
"Groundling" = "groundling\nsurvival",
"Fledgling" = "fledgling\nsurvival",
"Adult" = "adult\nsurvival")) +
scale_y_continuous(limits = c(-0.65, 0.65),
breaks = c(-0.5, -0.25, 0, 0.25, 0.5),
expand = c(0, 0), position = "left",
labels = c("-0.50","-0.25",
expression(phantom("-")*"0.00"),
expression(phantom("-")*"0.25"),
expression(phantom("-")*"0.50"))) +
xlab("Life stage") +
ylab("Sex bias\n")
Derived Adult Sex Ratio
ASR_boot <-
bind_rows(BC_hazard_rate_boot_tidy$ASR_ests_boot,
WBC_hazard_rate_boot_tidy$ASR_ests_boot) %>%
mutate(species = factor(species, levels = c("BC", "WBC")))
CI <- 0.95
ASR_boot_summary <-
ASR_boot %>%
dplyr::group_by(species) %>%
dplyr::summarise(ucl_ASR = stats::quantile(M_Adult, (1 - CI)/2, na.rm = TRUE),
lcl_ASR = stats::quantile(M_Adult, 1 - (1 - CI)/2, na.rm = TRUE),
avg_ASR = mean(M_Adult),
med_ASR = median(M_Adult),
max_ASR = max(M_Adult),
min_ASR = min(M_Adult))
Figure_2b <-
ggplot() +
annotate("rect", xmin=0, xmax=0.5, ymin=0, ymax=610, alpha=0.6,
fill= pull(ggthemes_data$wsj$palettes$colors6[3,2])) +
annotate("rect", xmin=0.5, xmax=1, ymin=0, ymax=610, alpha=0.6,
fill= pull(ggthemes_data$wsj$palettes$colors6[2,2])) +
annotate("text", x = c(0.05), y = c(305),
label = c("\u2640"), size = 4, colour = "grey10",
family="Menlo", vjust = c(0), hjust = c(0.5)) +
annotate("text", x = c(0.95), y = c(305),
label = c("\u2642"), size = 4, colour = "grey10",
family="Menlo", vjust = c(1), hjust = c(0.5)) +
geom_histogram(binwidth = 0.01, data = ASR_boot, aes(x = M_Adult), fill = "grey30") +
geom_errorbarh(data = ASR_boot_summary, aes(y = 600, x = lcl_ASR, xmin = lcl_ASR, xmax = ucl_ASR),
color = "black", size = 0.3, linetype = "solid") +
coord_flip() +
facet_grid(. ~ species) +
theme_bw() +
theme(text = element_text(family="Verdana"),
legend.position = "none",
axis.title.x = element_text(size=7, vjust=-0.1),
axis.text.x = element_text(size=6, angle = 45, hjust = 1, colour = "black"),
axis.title.y = element_text(size=7, hjust=0.5, vjust = -2.75, margin = margin(0, 11, 0, 0)),
axis.text.y = element_text(size=6, colour = "black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
axis.ticks.length = unit(0.1, "cm"),
axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
panel.border = element_blank(),
plot.margin = unit(c(0.2, 1.35, 0.405, 0.08), "cm"),
panel.spacing = unit(0.3, "lines"),
strip.background = element_blank(),
strip.text = element_blank()) +
ylab("Bootstrap frequency") +
xlab("Adult sex ratio\n(proportion \u2642, 95% CI)") +
scale_x_continuous(limits = c(0, 1), expand = c(0, 0)) +
scale_y_continuous(limits = c(0, 610), expand = c(0, 0), breaks=c(0, 100, 200, 300, 400, 500))
Life-table Response Experiment
Functions for the LTRE analysis
sensitivity_analysis()
Sensitivity analysis function takes the vital rate summary of the bootstrap procedure and conducts a perturbation analysis on each rate to assess how a proportional change in a given vital rate changes the ASR arguments: - vital_rate_summary: - matrix_str: - h: - k: - HSR: - niter: - ASR: - lambda:
sensitivity_analysis <-
function(vital_rate_summary, matrix_str, h = 1, k = 4,
HSR,
niter = 1000, ASR, lambda){
# make a list of all parameters
vr <-
list(F_Nestling_survival = vital_rate_summary$F_Nestling_survival,
F_Groundling_survival = vital_rate_summary$F_Groundling_survival,
F_Fledgling_survival = vital_rate_summary$F_Fledgling_survival,
F_Adult_survival = vital_rate_summary$F_Adult_survival,
M_Nestling_survival = vital_rate_summary$M_Nestling_survival,
M_Groundling_survival = vital_rate_summary$M_Groundling_survival,
M_Fledgling_survival = vital_rate_summary$M_Fledgling_survival,
M_Adult_survival = vital_rate_summary$M_Adult_survival)
# number of stages in the matrix
no_stages <- sqrt(length(matrix_str))
# Define plover life-stages of the Ceuta snowy plover matrix model
stages <- c("F_1st_year", "F_Adult", "M_1st_year", "M_Adult")
# an empty t by x matrix
stage <- matrix(numeric(no_stages * niter), nrow = no_stages)
# an empty t vector to store the population sizes
pop <- numeric(niter)
# dataframe to store the perturbation results
ASR_pert_results <-
data.frame(parameter = c("F_Nestling_survival", "F_Groundling_survival",
"F_Fledgling_survival", "F_Adult_survival",
"M_Nestling_survival", "F_Groundling_survival",
"M_Fledgling_survival", "M_Adult_survival",
"h", "k", "HSR", "ISR"),
sensitivities = numeric(length(vr) + 4),
elasticities = numeric(length(vr) + 4))
lambda_pert_results <-
data.frame(parameter = c("F_Nestling_survival", "F_Groundling_survival",
"F_Fledgling_survival", "F_Adult_survival",
"M_Nestling_survival", "F_Groundling_survival",
"M_Fledgling_survival", "M_Adult_survival",
"h", "k", "HSR", "ISR"),
sensitivities = numeric(length(vr) + 4),
elasticities = numeric(length(vr) + 4))
# specifiy how many survival rates there are
n <- length(vr)
# create vectors of perturbations to test on parameters of the matrix model
vr_nums <- seq(0, 1, 0.01) # proportional changes in survival and HSR (i.e., between 0 an 1)
h_nums <- seq(0, 2, 0.02) # proportional changes in h index (i.e., between 0 and 2)
k_nums <- seq(3, 5, 0.02) # proportional changes in k (i.e, between 3 and 5)
# create empty dataframes to store the perturbation results for ASR and lambda
vr_pert_ASR <- matrix(numeric(n * length(vr_nums)),
ncol = n, dimnames = list(vr_nums, names(vr)))
h_pert_ASR <- matrix(numeric(length(h_nums)),
ncol = 1, dimnames = list(h_nums, "h"))
k_pert_ASR <- matrix(numeric(length(k_nums)),
ncol = 1, dimnames = list(k_nums, "k"))
HSR_pert_ASR <- matrix(numeric(length(vr_nums)),
ncol = 1, dimnames = list(vr_nums, "HSR"))
vr_pert_lambda <- matrix(numeric(n * length(vr_nums)),
ncol = n, dimnames = list(vr_nums, names(vr)))
h_pert_lambda <- matrix(numeric(length(h_nums)),
ncol = 1, dimnames = list(h_nums, "h"))
k_pert_lambda <- matrix(numeric(length(k_nums)),
ncol = 1, dimnames = list(k_nums, "k"))
HSR_pert_lambda <- matrix(numeric(length(vr_nums)),
ncol = 1, dimnames = list(vr_nums, "HSR"))
#### perturbation of survival rates ####
for (g in 1:n) { # pick a column (i.e., a variable)
vr2 <- vr # reset the vital rates to the original
for (i in 1:length(vr_nums)) # pick a perturbation level
{
vr2[[g]] <- vr_nums[i] # specify the vital rate with the new perturbation level
A <- matrix(sapply(matrix_str, eval, vr2, NULL),
nrow = sqrt(length(matrix_str)), byrow=TRUE,
dimnames = list(stages, stages)) # build the matrix with the new value
# reset the starting stage distribution for simulation (all with 10 individuals)
m <- rep(10, no_stages)
for (j in 1:niter) { # project the matrix through t iteration
# stage distribution at time t
stage[,j] <- m
# population size at time t
pop[j] <- sum(m)
# number of male adults at time t
M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
# number of female adults at time t
F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
# Female freq-dep fecundity of Female chicks
A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - HSR) )
# Female freq-dep fecundity of Male chicks
A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * HSR)
# Male freq-dep fecundity of Female chicks
A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - HSR))
# Male freq-dep fecundity of Male chicks
A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * HSR)
# define the new n (i.e., new stage distribution at time t)
m <- A %*% (m)
}
# define rownames of stage matrix
rownames(stage) <- rownames(A)
# define colnames of stage matrix
colnames(stage) <- 0:(niter - 1)
# calculate the proportional stable stage distribution
stage <- apply(stage, 2, function(x) x / sum(x))
# define stable stage as the last stage
stable.stage <- stage[, niter]
# calc ASR as the proportion of the adult stable stage class that is male
vr_pert_ASR[i, g] <- stable.stage[no_stages] / (stable.stage[no_stages/2] +
stable.stage[no_stages])
# calc lambda as the pop change in the counts of the last two iterations
vr_pert_lambda[i, g] <- pop[niter]/pop[niter - 1]
}
# get the spline function of ASR
spl_ASR <-
smooth.spline(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g] ~
names(vr_pert_ASR[which(!is.na(vr_pert_ASR[, g])), g]))
# estimate the slope of the tangent of the spline at the vital rate
ASR_pert_results[g, 2] <- predict(spl_ASR, x = vr[[g]], deriv = 1)$y
# re-scale sensitivity into elasticity
ASR_pert_results[g, 3] <- vr[[g]] / ASR * ASR_pert_results[g, 2]
# do the same steps but for lambda
spl_lambda <-
smooth.spline(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g] ~
names(vr_pert_lambda[which(!is.na(vr_pert_lambda[, g])), g]))
lambda_pert_results[g, 2] <- predict(spl_lambda, x = vr[[g]], deriv = 1)$y
lambda_pert_results[g, 3] <- vr[[g]] / lambda * lambda_pert_results[g, 2]
}
#### perturbation of the h index parameter ####
for (i in 1:length(h_nums)) # pick a perturbation level
{
A <- matrix(sapply(matrix_str, eval, vr, NULL),
nrow = sqrt(length(matrix_str)), byrow=TRUE,
dimnames = list(stages, stages))
# reset the starting stage distribution for simulation (all with 10 individuals)
m <- rep(10, no_stages)
for (j in 1:niter) { # project the matrix through t iteration
# stage distribution at time t
stage[,j] <- m
# population size at time t
pop[j] <- sum(m)
# number of male adults at time t
M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
# number of female adults at time t
F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
# Female freq-dep fecundity of Female chicks
A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * (1 - HSR) )
# Female freq-dep fecundity of Male chicks
A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h_nums[i])) * HSR)
# Male freq-dep fecundity of Female chicks
A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * (1 - HSR))
# Male freq-dep fecundity of Male chicks
A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h_nums[i])) * HSR)
# define the new n (i.e., new stage distribution at time t)
m <- A %*% (m)
}
# define rownames of stage matrix
rownames(stage) <- rownames(A)
# define colnames of stage matrix
colnames(stage) <- 0:(niter - 1)
# calculate the proportional stable stage distribution
stage <- apply(stage, 2, function(x) x / sum(x))
# define stable stage as the last stage
stable.stage <- stage[, niter]
# calc ASR as the proportion of the adult stable stage class that is male
h_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages / 2] +
stable.stage[no_stages])
# calc lambda as the pop change in the counts of the last two iterations
h_pert_lambda[i, ] <- pop[niter]/pop[niter - 1]
}
# get the spline function of ASR
spl_ASR <-
smooth.spline(h_pert_ASR[which(!is.na(h_pert_ASR)), 1] ~
names(h_pert_ASR[which(!is.na(h_pert_ASR)), ]))
# estimate the slope of the tangent of the spline at the vital rate
ASR_pert_results[n + 1, 2] <- predict(spl_ASR, x = h, deriv = 1)$y
# re-scale sensitivity into elasticity
ASR_pert_results[n + 1, 3] <- h / ASR * ASR_pert_results[n + 1, 2]
# do the same steps but for lambda
spl_lambda <-
smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~
names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
lambda_pert_results[n + 1, 2] <- predict(spl_lambda, x = h, deriv = 1)$y
lambda_pert_results[n + 1, 3] <- h / lambda * lambda_pert_results[n + 1, 2]
#### perturbation of k parameter ####
for (i in 1:length(k_nums)) # pick a perturbation level
{
A <- matrix(sapply(matrix_str, eval, vr, NULL),
nrow = sqrt(length(matrix_str)), byrow=TRUE,
dimnames = list(stages, stages))
# reset the starting stage distribution for simulation (all with 10 individuals)
m <- rep(10, no_stages)
for (j in 1:niter) { # project the matrix through t iteration
# stage distribution at time t
stage[,j] <- m
# population size at time t
pop[j] <- sum(m)
# number of male adults at time t
M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
# number of female adults at time t
F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
# Female freq-dep fecundity of Female chicks
A["F_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * (1 - HSR) )
# Female freq-dep fecundity of Male chicks
A["M_1st_year", "F_Adult"] <- ((k_nums[i] * M2) / (M2 + (F2 / h)) * HSR)
# Male freq-dep fecundity of Female chicks
A["F_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * (1 - HSR))
# Male freq-dep fecundity of Male chicks
A["M_1st_year", "M_Adult"] <- ((k_nums[i] * F2) / (M2 + (F2 / h)) * HSR)
# define the new n (i.e., new stage distribution at time t)
m <- A %*% (m)
}
# define rownames of stage matrix
rownames(stage) <- rownames(A)
# define colnames of stage matrix
colnames(stage) <- 0:(niter - 1)
# calculate the proportional stable stage distribution
stage <- apply(stage, 2, function(x) x / sum(x))
# define stable stage as the last stage
stable.stage <- stage[, niter]
# calc ASR as the proportion of the adult stable stage class that is male
k_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] +
stable.stage[no_stages])
# calc lambda as the pop change in the counts of the last two iterations
k_pert_lambda[i, ] <- pop[niter]/pop[niter - 1]
}
# get the spline function of ASR
spl_ASR <-
smooth.spline(k_pert_ASR[which(!is.na(k_pert_ASR)), 1] ~
names(k_pert_ASR[which(!is.na(k_pert_ASR)), ]))
# estimate the slope of the tangent of the spline at the vital rate
ASR_pert_results[n + 2, 2] <- predict(spl_ASR, x = k, deriv = 1)$y
# re-scale sensitivity into elasticity
ASR_pert_results[n + 2, 3] <- k / ASR * ASR_pert_results[n + 2, 2]
# do the same steps but for lambda
spl_lambda <-
smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~
names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
lambda_pert_results[n + 2, 2] <- predict(spl_lambda, x = k, deriv = 1)$y
lambda_pert_results[n + 2, 3] <- k / lambda * lambda_pert_results[n + 2, 2]
#### perturbation of HSR ####
for (i in 1:length(vr_nums)) # pick a perturbation level
{
A <- matrix(sapply(matrix_str, eval, vr, NULL),
nrow = sqrt(length(matrix_str)), byrow=TRUE,
dimnames = list(stages, stages))
# reset the starting stage distribution for simulation (all with 10 individuals)
m <- rep(10, no_stages)
for (j in 1:niter) { # project the matrix through t iteration
# stage distribution at time t
stage[,j] <- m
# population size at time t
pop[j] <- sum(m)
# number of male adults at time t
M2 <- stage[4, j] * A["M_Adult", "M_Adult"]
# number of female adults at time t
F2 <- stage[2, j] * A["F_Adult", "F_Adult"]
# Female freq-dep fecundity of Female chicks
A["F_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * (1 - vr_nums[i]) )
# Female freq-dep fecundity of Male chicks
A["M_1st_year", "F_Adult"] <- ((k * M2) / (M2 + (F2 / h)) * vr_nums[i])
# Male freq-dep fecundity of Female chicks
A["F_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * (1 - vr_nums[i]))
# Male freq-dep fecundity of Male chicks
A["M_1st_year", "M_Adult"] <- ((k * F2) / (M2 + (F2 / h)) * vr_nums[i])
# define the new n (i.e., new stage distribution at time t)
m <- A %*% (m)
}
# define rownames of stage matrix
rownames(stage) <- rownames(A)
# define colnames of stage matrix
colnames(stage) <- 0:(niter - 1)
# calculate the proportional stable stage distribution
stage <- apply(stage, 2, function(x) x / sum(x))
# define stable stage as the last stage
stable.stage <- stage[, niter]
# calc ASR as the proportion of the adult stable stage class that is male
HSR_pert_ASR[i,] <- stable.stage[no_stages] / (stable.stage[no_stages/2] +
stable.stage[no_stages])
# calc lambda as the pop change in the counts of the last two iterations
HSR_pert_lambda[i, ] <- pop[niter] / pop[niter - 1]
}
# get the spline function of ASR
spl_ASR <-
smooth.spline(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), 1] ~
names(HSR_pert_ASR[which(!is.na(HSR_pert_ASR)), ]))
# estimate the slope of the tangent of the spline at the vital rate
ASR_pert_results[n + 3, 2] <- predict(spl_ASR, x = HSR, deriv = 1)$y
# re-scale sensitivity into elasticity
ASR_pert_results[n + 3, 3] <- HSR / ASR * ASR_pert_results[n + 3, 2]
# do the same steps but for lambda
spl_lambda <-
smooth.spline(h_pert_lambda[which(!is.na(h_pert_lambda)), 1] ~
names(h_pert_lambda[which(!is.na(h_pert_lambda)), ]))
lambda_pert_results[n + 3, 2] <- predict(spl_lambda, x = HSR, deriv = 1)$y
lambda_pert_results[n + 3, 3] <- HSR/lambda * lambda_pert_results[n + 3, 2]
#### store all results into a list ----
result <- list(ASR_pert_results = ASR_pert_results,
lambda_pert_results = lambda_pert_results,
HSR_pert_ASR = HSR_pert_ASR,
HSR_pert_lambda = HSR_pert_lambda,
k_pert_ASR = k_pert_ASR,
k_pert_lambda = k_pert_lambda,
h_pert_ASR = h_pert_ASR,
h_pert_lambda = h_pert_lambda,
vr_pert_ASR = vr_pert_ASR,
vr_pert_lambda = vr_pert_lambda)
}matrix_structure
matrix_structure is an expression that determines the structure of the two-sex matrix from a set of vital rates
matrix_structure <- expression(
# top row of matrix
0, NA, 0, NA,
# second row of matrix
(F_Nestling_survival * F_Groundling_survival * F_Fledgling_survival), F_Adult_survival, 0, 0,
# third row of matrix
0, NA, 0, NA,
# fourth row of matrix
0, 0, (M_Nestling_survival * M_Groundling_survival * M_Fledgling_survival), M_Adult_survival
)make_treat_matrix()
make_treat_matrix() fucntion makes a treatment matrix from a summary of the bootstrapped vital rates
arguments:
- survival_rates_boot_summary:
- species_name:
- h:
- k:
- HSR:
make_treat_matrix <-
function(survival_rates_boot_summary,
species_name,
h,
k,
HSR){
list(F_Nestling_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_nestling_survival")[, 4],
F_Groundling_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_groundling_survival")[, 4],
F_Fledgling_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_fledgling_survival")[, 4],
F_Adult_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_adult_survival")[, 4],
M_Nestling_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_nestling_survival")[, 4],
M_Groundling_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_groundling_survival")[, 4],
M_Fledgling_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_fledgling_survival")[, 4],
M_Adult_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_adult_survival")[, 4],
Egg_survival = filter(survival_rates_boot_summary,
species == species_name & vital_rate == "NA_egg_survival")[, 4],
# Define h (harem size, h = 1 is monogamy) and k (clutch size)
h = h,
k = k,
# Define primary sex ratio
HSR = HSR)
}make_mprime_matrix()
make_mprime_matrix() fucntion makes a prime-matrix (i.e., a matrix halfway between the treatment matrix and the unbiased matrix) from a set of vital rate summaries
arguments:
- survival_rates_boot_summary:
- species_name:
- h:
- k:
- HSR:
- sex:
make_mprime_matrix <-
function(survival_rates_boot_summary,
species_name,
h,
k,
HSR,
sex){
if(sex == "male"){
list(F_Nestling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
F_Groundling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
F_Fledgling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
F_Adult_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_adult_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
M_Nestling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_nestling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
M_Groundling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_groundling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
M_Fledgling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_fledgling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
M_Adult_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_adult_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
Egg_survival = filter(survival_rates_boot_summary,
species == species_name)[15, 4],
# Define h (harem size, h = 1 is monogamy) and k (clutch size)
h = (h + 1) / 2,
k = k,
# Define primary sex ratio
HSR = (HSR + 0.5) / 2)
}
else{
list(F_Nestling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_nestling_survival")[, 4]) / 2,
F_Groundling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_groundling_survival")[, 4]) / 2,
F_Fledgling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_fledgling_survival")[, 4]) / 2,
F_Adult_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_adult_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_adult_survival")[, 4]) / 2,
M_Nestling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_nestling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_nestling_survival")[, 4]) / 2,
M_Groundling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_groundling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_groundling_survival")[, 4]) / 2,
M_Fledgling_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_fledgling_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_fledgling_survival")[, 4]) / 2,
M_Adult_survival = (filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Female_adult_survival")[, 4] +
filter(survival_rates_boot_summary,
species == species_name & vital_rate == "Male_adult_survival")[, 4]) / 2,
Egg_survival = filter(survival_rates_boot_summary,
species == species_name)[15, 4],
# Define h (harem size, h = 1 is monogamy) and k (clutch size)
h = (h + 1) / 2,
k = k,
# Define primary sex ratio
HSR = (HSR + 0.5) / 2)
}
}LTRE_analysis()
LTRE_analysis() compares the relative difference in ASR response of an unbiased two-sex matrix vs. the observed sex-specific matrix
arguments:
- Mprime_sens:
- matrix_str:
- vital_rates:
- species_name:
- sex:
LTRE_analysis <-
function(Mprime_sens, matrix_str, vital_rates, species_name, sex){
# make empty dataframes to store LTRE results for ASR and lambda
LTRE_ASR <-
data.frame(parameter = c("Nestling survival", "Groundling survival",
"Fledgling survival", "Adult survival",
"Hatching sex ratio",
"Mating system"),
contribution = numeric(6))
LTRE_lambda <-
data.frame(parameter = c("Nestling survival", "Groundling survival",
"Fledgling survival", "Adult survival",
"Hatching sex ratio",
"Mating system"),
contribution = numeric(6))
# run a for loop to extract the parameter contributions
if (sex == "male"){
# male rates scenario
for(i in 1:nrow(LTRE_ASR))
{
LTRE_ASR[i, 2] <-
# survival rates
ifelse(i < 5, (vital_rates[[i + 4]] - vital_rates[[i]]) *
Mprime_sens$ASR_pert_results$sensitivities[i + 4],
# HSR
ifelse(i == 5, (vital_rates[[12]] - 0.5) *
Mprime_sens$ASR_pert_results$sensitivities[11],
# mating system
(1 - vital_rates[[10]]) *
Mprime_sens$ASR_pert_results$sensitivities[9]))
}
for(i in 1:nrow(LTRE_lambda))
{
# survival rates
ifelse(i < 5, (vital_rates[[i + 4]] - vital_rates[[i]]) *
Mprime_sens$lambda_pert_results$sensitivities[i + 4],
# HSR
ifelse(i == 5, (vital_rates[[12]] - 0.5) *
Mprime_sens$lambda_pert_results$sensitivities[11],
# mating system
(1 - vital_rates[[10]]) *
Mprime_sens$lambda_pert_results$sensitivities[9]))
}
}
else{
# female rates scenario
for(i in 1:nrow(LTRE_ASR))
{
LTRE_ASR[i, 2] <-
# survival rates
ifelse(i < 5, (vital_rates[[i]] - vital_rates[[i + 4]]) *
Mprime_sens$ASR_pert_results$sensitivities[i],
# HSR
ifelse(i == 5, (vital_rates[[12]] - 0.5) *
Mprime_sens$ASR_pert_results$sensitivities[11],
# mating system
(1 - vital_rates[[10]]) *
Mprime_sens$ASR_pert_results$sensitivities[9]))
}
for(i in 1:nrow(LTRE_lambda))
{
# survival rates
ifelse(i < 5, (vital_rates[[i]] - vital_rates[[i + 4]]) *
Mprime_sens$lambda_pert_results$sensitivities[i],
# HSR
ifelse(i == 5, (vital_rates[[12]] - 0.5) *
Mprime_sens$lambda_pert_results$sensitivities[11],
# mating system
(1 - vital_rates[[10]]) *
Mprime_sens$lambda_pert_results$sensitivities[9]))
}
}
LTRE_ASR$parameter <-
factor(LTRE_ASR$parameter, levels = c("Adult survival",
"Fledgling survival",
"Groundling survival",
"Nestling survival",
"Hatching sex ratio",
"Mating system"))
LTRE_lambda$parameter <-
factor(LTRE_lambda$parameter, levels = c("Adult survival",
"Fledgling survival",
"Groundling survival",
"Nestling survival",
"Hatching sex ratio",
"Mating system"))
LTRE_ASR$model <- "ASR"
LTRE_ASR$species <- species_name
LTRE_lambda$model <- "lambda"
LTRE_lambda$species <- species_name
LTRE_results <- list(LTRE_ASR = LTRE_ASR,
LTRE_lambda = LTRE_lambda)
LTRE_results
}LTRE_contributions_check()
LTRE_contributions_check() takes the results of the LTRE analysis and checks the total contribution that each vital rate has on ASR. In general, the sum of the LTRE contribution should roughly equal the difference in ASR between the observed scenario and the unbiased scenario.
arguments:
- M_matrix_vital_rates:
- Mprime_sensitivities:
- M_matrix_ASR:
- scenario:
LTRE_contributions_check <-
function(M_matrix_vital_rates, Mprime_sensitivities, M_matrix_ASR, scenario){
if (scenario == "male")
contribution_sum <-
sum(
(M_matrix_vital_rates[[1]] - M_matrix_vital_rates[[5]]) *
Mprime_sensitivities[1],
(M_matrix_vital_rates[[2]] - M_matrix_vital_rates[[6]]) *
Mprime_sensitivities[2],
(M_matrix_vital_rates[[3]] - M_matrix_vital_rates[[7]]) *
Mprime_sensitivities[3],
(M_matrix_vital_rates[[4]] - M_matrix_vital_rates[[8]]) *
Mprime_sensitivities[4],
(M_matrix_vital_rates[[5]] - M_matrix_vital_rates[[5]]) *
Mprime_sensitivities[1],
(M_matrix_vital_rates[[6]] - M_matrix_vital_rates[[6]]) *
Mprime_sensitivities[2],
(M_matrix_vital_rates[[7]] - M_matrix_vital_rates[[7]]) *
Mprime_sensitivities[3],
(M_matrix_vital_rates[[8]] - M_matrix_vital_rates[[8]]) *
Mprime_sensitivities[4],
(1 - M_matrix_vital_rates[[10]]) *
Mprime_sensitivities[9],
(M_matrix_vital_rates[[12]] - 0.5) * Mprime_sensitivities[11],
(M_matrix_vital_rates[[13]] - 0.5) * Mprime_sensitivities[12]
)
else
contribution_sum <-
sum(
(M_matrix_vital_rates[[5]] - M_matrix_vital_rates[[1]]) *
Mprime_sensitivities[5],
(M_matrix_vital_rates[[6]] - M_matrix_vital_rates[[2]]) *
Mprime_sensitivities[6],
(M_matrix_vital_rates[[7]] - M_matrix_vital_rates[[3]]) *
Mprime_sensitivities[7],
(M_matrix_vital_rates[[8]] - M_matrix_vital_rates[[4]]) *
Mprime_sensitivities[8],
(M_matrix_vital_rates[[1]] - M_matrix_vital_rates[[1]]) *
Mprime_sensitivities[5],
(M_matrix_vital_rates[[2]] - M_matrix_vital_rates[[2]]) *
Mprime_sensitivities[6],
(M_matrix_vital_rates[[3]] - M_matrix_vital_rates[[3]]) *
Mprime_sensitivities[7],
(M_matrix_vital_rates[[4]] - M_matrix_vital_rates[[4]]) *
Mprime_sensitivities[8],
(1 - M_matrix_vital_rates[[10]]) *
Mprime_sensitivities[9],
(M_matrix_vital_rates[[12]] - 0.5) * Mprime_sensitivities[11],
(M_matrix_vital_rates[[13]] - 0.5) * Mprime_sensitivities[12]
)
ASR_bias <- abs(M_matrix_ASR - 0.5)
absolute_difference <- abs(ASR_bias) - abs(contribution_sum)
return(list(contribution_sum = as.vector(contribution_sum),
ASR_bias = as.vector(ASR_bias),
absolute_difference = as.vector(absolute_difference)))
}LTRE: Analysis
BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <-
as.factor(BC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)
WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter <-
as.factor(WBC_hazard_rate_boot_tidy$vital_rate_ests_boot$iter)
# summarize the vital rates
survival_rates_boot_summary <-
bind_rows(BC_hazard_rate_boot_tidy$vital_rate_ests_boot, WBC_hazard_rate_boot_tidy$vital_rate_ests_boot) %>%
mutate(vital_rate = paste(sex, stage, rate, sep = "_")) %>%
Rmisc::summarySE(.,
measurevar = "value",
groupvars = c("vital_rate", "species"),
conf.interval = 0.95) %>%
arrange(species, vital_rate)
BC_VR_treat <-
make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
species_name = "BC")
BC_VR_mprime_male <-
make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
species_name = "BC",
sex = "male")
BC_VR_mprime_female <-
make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
h = 1/pull(filter(parameter_distributions, species == "BC" & trait == "mating_system"), mean),
HSR = pull(filter(parameter_distributions, species == "BC" & trait == "hatching_sex_ratio"), mean),
k = pull(filter(parameter_distributions, species == "BC" & trait == "clutch_size"), mean),
species_name = "BC",
sex = "female")
WBC_VR_treat <-
make_treat_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
species_name = "WBC")
WBC_VR_mprime_male <-
make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
species_name = "WBC",
sex = "male")
WBC_VR_mprime_female <-
make_mprime_matrix(survival_rates_boot_summary = survival_rates_boot_summary,
h = 1/pull(filter(parameter_distributions, species == "WBC" & trait == "mating_system"), mean),
HSR = pull(filter(parameter_distributions, species == "WBC" & trait == "hatching_sex_ratio"), mean),
k = pull(filter(parameter_distributions, species == "WBC" & trait == "clutch_size"), mean),
species_name = "WBC",
sex = "female")
BC_treatment_matrix <- coucal_matrix(BC_VR_treat)
BC_M_prime_matrix_male <- coucal_matrix(BC_VR_mprime_male)
BC_M_prime_matrix_female <- coucal_matrix(BC_VR_mprime_female)
WBC_treatment_matrix <- coucal_matrix(WBC_VR_treat)
WBC_M_prime_matrix_male <- coucal_matrix(WBC_VR_mprime_male)
WBC_M_prime_matrix_female <- coucal_matrix(WBC_VR_mprime_female)
BC_treatment_ASR_analysis <-
matrix_ASR(M = BC_treatment_matrix,
h = BC_VR_treat$h,
k = BC_VR_treat$k,
HSR = BC_VR_treat$HSR,
iterations = 100,
num_boot = 1,
iter_add = 1,
species = "BC")
BC_ASR_treat <- BC_treatment_ASR_analysis$ASR
BC_ASR_treat M_Adult
0.6901214 WBC_treatment_ASR_analysis <-
matrix_ASR(M = WBC_treatment_matrix,
h = WBC_VR_treat$h,
k = WBC_VR_treat$k,
HSR = WBC_VR_treat$HSR,
iterations = 100,
num_boot = 1,
iter_add = 1,
species = "BC")
WBC_ASR_treat <- WBC_treatment_ASR_analysis$ASR
WBC_ASR_treat M_Adult
0.4936118 BC_M_prime_ASR_analysis_male <-
matrix_ASR(M = BC_M_prime_matrix_male,
h = BC_VR_mprime_male$h,
k = BC_VR_mprime_male$k,
HSR = BC_VR_mprime_male$HSR,
iterations = 100,
num_boot = 1,
iter_add = 1,
species = "BC")
BC_ASR_mprime_male <- BC_M_prime_ASR_analysis_male$ASR
BC_ASR_mprime_male M_Adult
0.5855401 BC_M_prime_ASR_analysis_female <-
matrix_ASR(M = BC_M_prime_matrix_female,
h = BC_VR_mprime_female$h,
k = BC_VR_mprime_female$k,
HSR = BC_VR_mprime_female$HSR,
iterations = 100,
num_boot = 1,
iter_add = 1,
species = "BC")
BC_ASR_mprime_female <- BC_M_prime_ASR_analysis_female$ASR
BC_ASR_mprime_female M_Adult
0.5797847 WBC_M_prime_ASR_analysis_male <-
matrix_ASR(M = WBC_M_prime_matrix_male,
h = WBC_VR_mprime_male$h,
k = WBC_VR_mprime_male$k,
HSR = WBC_VR_mprime_male$HSR,
iterations = 100,
num_boot = 1,
iter_add = 1,
species = "BC")
WBC_ASR_mprime_male <- WBC_M_prime_ASR_analysis_male$ASR
WBC_ASR_mprime_male M_Adult
0.4951691 WBC_M_prime_ASR_analysis_female <-
matrix_ASR(M = WBC_M_prime_matrix_female,
h = WBC_VR_mprime_female$h,
k = WBC_VR_mprime_female$k,
HSR = WBC_VR_mprime_female$HSR,
iterations = 100,
num_boot = 1,
iter_add = 1,
species = "BC")
WBC_ASR_mprime_female <- WBC_M_prime_ASR_analysis_female$ASR
WBC_ASR_mprime_female M_Adult
0.4983036 BC_lambda_treat <- BC_treatment_ASR_analysis$lambda
BC_lambda_treat[1] 0.8686866BC_lambda_mprime_male <- BC_M_prime_ASR_analysis_male$lambda
BC_lambda_mprime_male[1] 1.008295BC_lambda_mprime_female <- BC_M_prime_ASR_analysis_female$lambda
BC_lambda_mprime_female[1] 0.8823659WBC_lambda_treat <- WBC_treatment_ASR_analysis$lambda
WBC_lambda_treat[1] 0.9134844WBC_lambda_mprime_male <- WBC_M_prime_ASR_analysis_male$lambda
WBC_lambda_mprime_male[1] 0.920894WBC_lambda_mprime_female <- WBC_M_prime_ASR_analysis_female$lambda
WBC_lambda_mprime_female[1] 0.9242015BC_treat_sensitivity_analysis <-
sensitivity_analysis(vital_rate_summary = BC_VR_treat,
matrix_str = matrix_structure,
h = BC_VR_treat$h,
k = BC_VR_treat$k,
HSR = BC_VR_treat$HSR,
niter = 100,
ASR = BC_ASR_treat,
lambda = BC_lambda_treat)
BC_Mprime_sensitivity_analysis_male <-
sensitivity_analysis(vital_rate_summary = BC_VR_mprime_male,
matrix_str = matrix_structure,
h = BC_VR_mprime_male$h,
k = BC_VR_mprime_male$k,
HSR = BC_VR_mprime_male$HSR,
niter = 100,
ASR = BC_ASR_mprime_male,
lambda = BC_lambda_mprime_male)
BC_Mprime_sensitivity_analysis_female <-
sensitivity_analysis(vital_rate_summary = BC_VR_mprime_female,
matrix_str = matrix_structure,
h = BC_VR_mprime_female$h,
k = BC_VR_mprime_female$k,
HSR = BC_VR_mprime_female$HSR,
niter = 100,
ASR = BC_ASR_mprime_female,
lambda = BC_lambda_mprime_female)
WBC_treat_sensitivity_analysis <-
sensitivity_analysis(vital_rate_summary = WBC_VR_treat,
matrix_str = matrix_structure,
h = WBC_VR_treat$h,
k = WBC_VR_treat$k,
HSR = WBC_VR_treat$HSR,
niter = 100,
ASR = WBC_ASR_treat,
lambda = WBC_lambda_treat)
WBC_Mprime_sensitivity_analysis_male <-
sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_male,
matrix_str = matrix_structure,
h = WBC_VR_mprime_male$h,
k = WBC_VR_mprime_male$k,
HSR = WBC_VR_mprime_male$HSR,
niter = 100,
ASR = WBC_ASR_mprime_male,
lambda = WBC_lambda_mprime_male)
WBC_Mprime_sensitivity_analysis_female <-
sensitivity_analysis(vital_rate_summary = WBC_VR_mprime_female,
matrix_str = matrix_structure,
h = WBC_VR_mprime_female$h,
k = WBC_VR_mprime_female$k,
HSR = WBC_VR_mprime_female$HSR,
niter = 100,
ASR = WBC_ASR_mprime_female,
lambda = WBC_lambda_mprime_female)BC_LTRE_male <-
LTRE_analysis(Mprime_sens = BC_Mprime_sensitivity_analysis_male,
matrix_str = matrix_str,
vital_rates = BC_VR_treat,
species_name = "BC",
sex = "male")
BC_LTRE_female <-
LTRE_analysis(Mprime_sens = BC_Mprime_sensitivity_analysis_female,
matrix_str = matrix_str,
vital_rates = BC_VR_treat,
species_name = "BC",
sex = "female")
WBC_LTRE_male <-
LTRE_analysis(Mprime_sens = WBC_Mprime_sensitivity_analysis_male,
matrix_str = matrix_str,
vital_rates = WBC_VR_treat,
species_name = "WBC",
sex = "male")
WBC_LTRE_female <-
LTRE_analysis(Mprime_sens = WBC_Mprime_sensitivity_analysis_female,
matrix_str = matrix_str,
vital_rates = WBC_VR_treat,
species_name = "WBC",
sex = "female")
LTRE_coucal_male_ASR <-
bind_rows(BC_LTRE_male$LTRE_ASR, WBC_LTRE_male$LTRE_ASR) %>%
mutate(sex = "male")
LTRE_coucal_female_ASR <-
bind_rows(BC_LTRE_female$LTRE_ASR, WBC_LTRE_female$LTRE_ASR) %>%
mutate(sex = "female")
LTRE_coucal_ASR <- rbind(LTRE_coucal_male_ASR, LTRE_coucal_female_ASR)
LTRE_coucal_ASR$parameter <-
factor(LTRE_coucal_ASR$parameter,
levels = c("Hatching sex ratio",
"Nestling survival",
"Groundling survival",
"Fledgling survival",
"Adult survival",
"Mating system"))LTRE Results
# custom color palette for the plotting of Juvenile and Adult stats
cbPalette <- c("#D9D9D9", "#D9D9D9", "#D9D9D9",
"#D9D9D9", "#A6A6A6", "#A6A6A6",
"#A6A6A6")
species_names <-
c('BC' = "Black Coucal",
'WBC' = "White-browed Coucal")
analysis_names <- c(
'male' = "Male Mo scenario",
'female' = "Female Mo scenario"
)
# plot the comparative LTRE results
vital_rate_theme <-
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.length = unit(0.1, "cm"),
panel.border = element_blank(),
panel.spacing.x = unit(0.3, "lines"),
panel.spacing.y = unit(0.7, "lines"),
strip.background = element_blank()
)
# theme_set(vital_rate_theme)
LTRE_ASR_plot_background <-
ggplot2::ggplot(data = LTRE_coucal_ASR,
aes(x = parameter, y = contribution, fill = parameter)) +
annotate("rect", xmin=-Inf, xmax=Inf, ymin=-0.15, ymax=0, alpha=0.6,
fill= sex_pal2[1]) +
annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=0.15, alpha=0.6,
fill= sex_pal2[2]) +
annotate("text", x = c(2), y = c(-0.14),
label = c("\u2640"), size = 4, family = "Menlo",
vjust = c(0), hjust = c(0.5), colour = "grey10") +
annotate("text", x = c(2), y = c(0.14),
label = c("\u2642"), size = 4, family = "Menlo",
vjust = c(1), hjust = c(0.5), colour = "grey10") +
facet_grid(sex ~ species,
labeller = labeller(.cols = species_names, .rows = analysis_names)) +
vital_rate_theme +
theme(axis.title.x = element_text(size = 7, vjust = -0.1, colour = "white"),
axis.text.x = element_text(size = 7, angle = 45, hjust = 1, colour = "white"),
axis.ticks.x = element_line(size = 0.2, colour = "white"),
strip.text.x = element_text(size = 6, colour = "white"),
axis.title.y = element_text(size = 7, hjust=0.5, vjust = 3.5, colour = "white"),
axis.text.y = element_text(size = 6, colour = "white"),
axis.ticks.y = element_blank(),
strip.text.y = element_text(size = 6, colour = "white"),
plot.margin = unit(c(0.2, 0.85,
0.78, 0.6), "cm")) +
scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
"Nestling survival" = expression(italic(S["n"])),
"Groundling survival" = expression(italic(S["g"])),
"Fledgling survival" = expression(italic(S["f"])),
"Adult survival" = expression(italic(phi["a"])),
"Mating system" = expression(italic("h")))) +
scale_y_continuous(limits = c(-0.15 , 0.15), expand = c(0, 0),
breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
ylab("Absolute contribution to adult sex ratio bias") +
xlab("Sex bias in parameter")
LTRE_ASR_plot <-
ggplot2::ggplot() +
geom_bar(data = LTRE_coucal_ASR,
aes(x = parameter, y = contribution, fill = parameter),
color = "black", stat = "identity", size = 0.2) +
facet_grid(sex ~ species,
labeller = labeller(.cols = species_names, .rows = analysis_names)) +
vital_rate_theme +
theme(axis.title.x = element_text(size=7, vjust=-0.1),
axis.text.x = element_text(size=7, angle = 45, hjust = 1),
axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
strip.text = element_text(size = 6, colour = "black", face = "italic"),
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
axis.title.y = element_text(size=7, hjust=0.5, vjust = 3.5),
axis.text.y = element_text(size=6),
axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
strip.text.y = element_text(size = 6),
plot.margin = unit(c(0.2, 0.85, 0.2, 0.6), "cm")) +
scale_fill_manual(values = cbPalette) +
scale_y_continuous(limits = c(-0.15 ,0.15), expand = c(0, 0),
breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
"Nestling survival" = expression(italic(S["n"])),
"Groundling survival" = expression(italic(S["g"])),
"Fledgling survival" = expression(italic(S["f"])),
"Adult survival" = expression(italic(phi["a"])),
"Mating system" = expression(italic("h")))) +
ylab("Absolute contribution to adult sex ratio bias") +
xlab("Sex bias in parameter")
# jpeg(filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_ASR_plot.jpeg",
# # compression = "none",
# width = 4.75,
# height = 3,
# units = "in",
# res = 600)
grid::grid.newpage()
pushViewport( viewport( layout = grid.layout( 1 , 1 , widths = unit( 1 , "npc" ) ) ) )
print(LTRE_ASR_plot_background, newpage = FALSE)
print(LTRE_ASR_plot, newpage = FALSE)
grid::popViewport()
# dev.off()# plot the comparative LTRE results
vital_rate_theme <-
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.position = "none",
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.length = unit(0.1, "cm"),
panel.border = element_blank(),
panel.spacing.x = unit(0.3, "lines"),
panel.spacing.y = unit(0.7, "lines"),
strip.background = element_blank()
)
LTRE_ASR_plot_background_female <-
ggplot2::ggplot(data = filter(LTRE_coucal_ASR, sex == "female"),
aes(x = parameter, y = contribution, fill = parameter)) +
annotate("rect", xmin=-Inf, xmax=Inf, ymin=-0.15, ymax=0, alpha=0.6,
fill= sex_pal2[1]) +
annotate("rect", xmin=-Inf, xmax=Inf, ymin=0, ymax=0.15, alpha=0.6,
fill= sex_pal2[2]) +
annotate("text", x = c(2), y = c(-0.14),
label = c("\u2640"), size = 4, family = "Menlo",
vjust = c(0), hjust = c(0.5), colour = "grey10") +
annotate("text", x = c(2), y = c(0.14),
label = c("\u2642"), size = 4, family = "Menlo",
vjust = c(1), hjust = c(0.5), colour = "grey10") +
facet_grid(. ~ species,
labeller = labeller(.cols = species_names)) +
vital_rate_theme +
theme(axis.title.x = element_text(size = 7, vjust = -0.1, colour = "white"),
axis.text.x = element_text(size = 7, angle = 45, hjust = 1, colour = "white"),
axis.ticks.x = element_line(size = 0.2, colour = "white"),
strip.text.x = element_text(size = 6, colour = "white", face = "italic"),
axis.title.y = element_text(size = 7, hjust=0.5, vjust = 3.5, colour = "white"),
axis.text.y = element_text(size = 6, colour = "white"),
axis.ticks.y = element_blank(),
strip.text.y = element_text(size = 6, colour = "white"),
plot.margin = unit(c(0.2, 0.85,
0.78, 0.6), "cm")) +
scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
"Nestling survival" = expression(italic(S["n"])),
"Groundling survival" = expression(italic(S["g"])),
"Fledgling survival" = expression(italic(S["f"])),
"Adult survival" = expression(italic(phi["a"])),
"Mating system" = expression(italic("h")))) +
scale_y_continuous(limits = c(-0.15 , 0.15), expand = c(0, 0),
breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
ylab("Absolute contribution to\nadult sex ratio bias") +
xlab("Sex bias in parameter")
LTRE_ASR_plot_female <-
ggplot2::ggplot() +
geom_bar(data = filter(LTRE_coucal_ASR, sex == "female"),
aes(x = parameter, y = contribution, fill = parameter),
color = "black", stat = "identity", size = 0.2) +
facet_grid(. ~ species,
labeller = labeller(.cols = species_names)) +
vital_rate_theme +
theme(axis.title.x = element_text(size=7, vjust=-0.1),
axis.text.x = element_text(size=7, angle = 45, hjust = 1, colour = "white"),
axis.ticks.x = element_line(size = 0.2, colour = "grey40"),
strip.text = element_text(size = 6, colour = "black", face = "italic"),
panel.background = element_rect(fill = "transparent", colour = NA),
plot.background = element_rect(fill = "transparent", colour = NA),
axis.title.y = element_text(size=7, hjust=0.5, vjust = 3.5),
axis.text.y = element_text(size=6),
axis.ticks.y = element_line(size = 0.2, colour = "grey40"),
strip.text.y = element_text(size = 6),
plot.margin = unit(c(0.2, 0.85, 0.2, 0.6), "cm")) +
scale_fill_manual(values = cbPalette) +
scale_y_continuous(limits = c(-0.15 ,0.15), expand = c(0, 0),
breaks = c(-0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15)) +
scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
"Nestling survival" = expression(italic(S["n"])),
"Groundling survival" = expression(italic(S["g"])),
"Fledgling survival" = expression(italic(S["f"])),
"Adult survival" = expression(italic(phi["a"])),
"Mating system" = expression(italic("h")))) +
ylab("Absolute contribution to\nadult sex ratio bias") +
xlab("Sex bias in parameter")
# jpeg(filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_ASR_plot_female.jpeg",
# # compression = "none",
# width = 4,
# height = 1.5,
# units = "in",
# res = 600)
grid.newpage()
pushViewport( viewport( layout = grid.layout( 1 , 1 , widths = unit( 1 , "npc" ) ) ) )
print(LTRE_ASR_plot_background_female, newpage = FALSE)
print(LTRE_ASR_plot_female, newpage = FALSE)
grid::popViewport()
# dev.off()# Determine how much larger the contribution of each vital rates is compared to juvenile survival
# groundling vs nestling:
LTRE_coucal_ASR_prop <-
filter(LTRE_coucal_ASR) %>%
dplyr::select(-model) %>%
mutate(parameter2 = as.character(parameter)) %>%
dplyr::rename(parameter1 = parameter)
prop_contribution_table <-
LTRE_coucal_ASR_prop %>%
tidyr::expand(., parameter1 = parameter1, parameter2 = parameter1) %>%
left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter1, contribution, sex, species),
by = c("parameter1")) %>%
dplyr::rename(contribution1 = contribution) %>%
left_join(., dplyr::select(LTRE_coucal_ASR_prop, parameter2, contribution, sex, species),
by = c("parameter2", "species", "sex")) %>%
dplyr::rename(contribution2 = contribution) %>%
mutate(prop_diff = abs(contribution1) / abs(contribution2),
prop_diff2 = ifelse(contribution1 < 0 & contribution2 < 0, -(contribution1/contribution2),
ifelse(contribution1 < 0 & contribution2 > 0, contribution1/contribution2,
ifelse(contribution1 > 0 & contribution2 < 0, -(contribution1/contribution2), contribution1/contribution2))),
parameter1 = factor(parameter1,
levels = c("Hatching sex ratio",
"Nestling survival",
"Groundling survival",
"Fledgling survival",
"Adult survival",
"Mating system")),
parameter2 = factor(parameter2,
levels = c("Hatching sex ratio",
"Nestling survival",
"Groundling survival",
"Fledgling survival",
"Adult survival",
"Mating system")))
LTRE_heatmap <-
prop_contribution_table %>%
filter(parameter1 %in% c("Mating system",
"Adult survival",
"Fledgling survival",
"Groundling survival",
"Nestling survival",
"Hatching sex ratio") &
parameter2 %in% c("Mating system",
"Adult survival",
"Fledgling survival",
"Groundling survival",
"Nestling survival",
"Hatching sex ratio")) %>%
ggplot(.,
aes(parameter1, parameter2, fill = prop_diff2)) +
geom_tile() +
facet_grid(sex ~ species,
labeller = labeller(.cols = species_names, .rows = analysis_names)) +
theme(
axis.text.x = element_text(size = 7, angle = 45, hjust = 1),
axis.text.y = element_text(size = 7),
axis.title.x = element_text(size = 7),
axis.title.y = element_text(size = 7),
legend.position = "bottom",
legend.box = "horizontal",
legend.text = element_text(size = 7),
legend.title = element_text(size = 7),
axis.ticks = element_blank(),
panel.border = element_rect(color = "grey40", fill = NA),
) +
scale_x_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
"Nestling survival" = expression(italic(S["n"])),
"Groundling survival" = expression(italic(S["g"])),
"Fledgling survival" = expression(italic(S["f"])),
"Adult survival" = expression(italic(phi["a"])),
"Mating system" = expression(italic("h")))) +
scale_y_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
"Nestling survival" = expression(italic(S["n"])),
"Groundling survival" = expression(italic(S["g"])),
"Fledgling survival" = expression(italic(S["f"])),
"Adult survival" = expression(italic(phi["a"])),
"Mating system" = expression(italic("h")))) +
scale_fill_gradient2(low = sex_pal2[1],
mid = "white",
high = sex_pal2[2],
name = "Proportional contribution of A compared to B") +#,
xlab("Parameter A") +
ylab("Parameter B")
LTRE_heatmap
# ggsave(LTRE_heatmap,
# filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_heatmap.jpeg",
# # compression = "none",
# width = 4,
# height = 4,
# units = "in",
# dpi = 600)
# LTRE_heatmap_female <-
# prop_contribution_table %>%
# mutate(prop_diff2 = ifelse(abs(prop_diff2) == 1, 0, prop_diff2)) %>%
# filter(sex == "female") %>%
# filter(parameter1 %in% c("Mating system",
# "Adult survival",
# "Fledgling survival",
# "Groundling survival",
# "Nestling survival",
# "Hatching sex ratio") &
# parameter2 %in% c("Mating system",
# "Adult survival",
# "Fledgling survival",
# "Groundling survival",
# "Nestling survival",
# "Hatching sex ratio")) %>%
# ggplot(.,
# aes(parameter1, parameter2, fill = prop_diff2)) +
# geom_tile() +
# facet_grid(. ~ species,
# labeller = labeller(.cols = species_names)) +
# theme(
# axis.text.x = element_text(size = 7, angle = 45, hjust = 1),
# axis.text.y = element_text(size = 7),
# axis.title.x = element_text(size = 7),
# axis.title.y = element_text(size = 7),
# legend.position = "bottom",
# legend.box = "horizontal",
# legend.text = element_text(size = 5),
# legend.title = element_text(size = 6),
# axis.ticks = element_line(size = 0.2, colour = "grey40"),
# panel.border = element_rect(color = "grey40", fill = NA),
# plot.margin = unit(c(0.2, 0.85, 0.2,
# 0.55), "cm"),
# strip.text = element_blank()
# ) +
# scale_x_discrete(labels = c("HSR" = "hatching\nsex ratio",
# "Nestling" = "nestling\nsurvival",
# "Groundling" = "groundling\nsurvival",
# "Fledgling" = "fledgling\nsurvival",
# "Adult" = "adult\nsurvival",
# "Mating system" = "Mating system")) +
# scale_y_discrete(labels = c("Hatching sex ratio" = expression(italic("\u03C1")),
# "Nestling survival" = expression(italic(S["n"])),
# "Groundling survival" = expression(italic(S["g"])),
# "Fledgling survival" = expression(italic(S["f"])),
# "Adult survival" = expression(italic(phi["a"])),
# "Mating system" = expression(italic("h")))) +
# scale_fill_gradient2(low = sex_pal2[1],
# mid = "white",
# high = sex_pal2[2],
# name = "Proportional contribution of A compared to B") +#,
# xlab("Sex-bias in parameter A") +
# ylab("Sex-bias in\nparameter B\n ")
# LTRE_heatmap_female
# ggsave(LTRE_heatmap_female,
# filename = "/Users/luketheduke2/ownCloud/coucal_ASR/R_project/products/figures/LTRE_heatmap_female.jpeg",
# # compression = "none",
# width = 4,
# height = 2.5,
# units = "in",
# dpi = 600)Supplementary Analysis 1: sex-specific juvenile movement
coucal_wp <-
read.csv("data/raw/juvies_2013_2020_waypoints.csv",
stringsAsFactors = FALSE, header = TRUE) %>%
# clean up a few inconsistencies in the data
mutate(timestamp = ifelse(ring_ID == "GN 76660" & date_dec == "20140502", "2014-05-02 10:09:00.000", timestamp),
location.long = ifelse(location.long == 34.68357, 34.08357, location.long)) %>%
filter(ring_ID != "GN 76726") %>%
filter(ring_ID != "GN 76728") %>%
filter(ring_ID != "GN 76732") %>%
filter(ring_ID != "GN 76733") %>%
filter(ring_ID != "GN 76777") %>% # both male and female??
filter(movebank_event.id != "16897537877") %>%
filter(movebank_event.id != "16897535798") %>%
filter(movebank_event.id != "16897537758") %>%
filter(movebank_event.id != "16897538205") %>%
# specify the timezone of the time stamp
mutate(timestamp = as.POSIXct(timestamp,
tz = "Africa/Nairobi"),
# make the serial number a 4 digit string
serial_number = as.character(str_pad(serial_number, 4, pad = "0"))) %>%
# make a julian date time variable for plotting
mutate(julian_date = as.numeric(format(timestamp, "%j"))) %>%
# remove observations older than 70 days of age
filter(age_days < 71) %>%
# consolidate columns of interest
dplyr::select(species, year, nest.nr, ring_ID, sex, timestamp, julian_date,
age_days, location.long, location.lat, species, site, location,
comment)
# specify coordinate columns
coordinates(coucal_wp) <- c("location.long","location.lat")
# define spatial projection as WGS84
proj4string(coucal_wp) <- CRS("+init=epsg:4326")# specify mapview options for viewing
mapviewOptions(basemaps = c("Esri.WorldImagery"),
layers.control.pos = "topright",
legend.pos = "bottomleft")
# check out the spatial data
mapview(coucal_wp[coucal_wp$ring_ID %in% "GN 76751",], zcol = "julian_date",
layer.name = "WBC radio tagging",
layers.control.pos = "topright")# make a spatial trajectory object of each animal's encounter history
tr_coucal_wp <-
as.ltraj(xy = sp::coordinates(coucal_wp),
date = coucal_wp$timestamp,
id = coucal_wp$ring_ID)
# convert the trajectory object back to a dataframe, rename columns, and
# calculate the cumulative distance traveled over the observation
tr_coucal_wp <-
ld(tr_coucal_wp) %>%
dplyr::rename(distance = dist,
dispersion = R2n,
cardinal_angle = abs.angle,
relative_angle = rel.angle,
ring_ID = id) %>%
group_by(ring_ID) %>%
mutate(cum_distance = cumsum(distance))
# extract the temporal summary info for each bird (min age, min date, max age)
# and merge to the trajectory dataframe
coucal_wp_clean <-
as.data.frame(coucal_wp) %>%
group_by(species, ring_ID, sex, nest.nr) %>%
dplyr::summarise(min_age = min(age_days),
min_date = min(timestamp),
max_age = max(age_days)) %>%
arrange(desc(max_age)) %>%
left_join(., tr_coucal_wp, by = "ring_ID") %>%
# calcuate the differences in days between the observation and minimum date
# to get time since first observation, convert from seconds to days, then add
# to the minimum age value to extract the age at a given observation
mutate(date_diff = date - min_date,
year = year(date)) %>%
mutate(date_diff = round(((as.numeric(date_diff)/60)/60)/24)) %>%
mutate(age = min_age + date_diff) %>%
ungroup() %>%
# remove all NA rows
filter(!is.na(cum_distance)) %>%
# calculate the temporal lag between observations
group_by(ring_ID) %>%
mutate(obs_diff = date_diff - lag(date_diff),
origin_x = x[which.min(date)],
origin_y = y[which.min(date)])
# for each bird calculate the daily distance from origin
coucal_wp_clean$distance_origin <-
sapply(1:nrow(coucal_wp_clean), function(i)
spDistsN1(as.matrix(coucal_wp_clean[i, c("x", "y")]),
as.matrix(coucal_wp_clean[i, c("origin_x", "origin_y")]),
longlat = TRUE))
# transform the distance metrics to meters
coucal_wp_clean <-
coucal_wp_clean %>%
mutate(distance = distance * 100000,
distance_origin = distance_origin * 1000,
cum_distance = cum_distance * 100000)
# check the distribution of distances
hist(log(coucal_wp_clean$distance))
range(coucal_wp_clean$distance)
boxplot.stats(coucal_wp_clean$distance)# extract individuals that mainly had observations spaced 2 days apart
two_day_rings <-
coucal_wp_clean %>%
dplyr::select(species, ring_ID, sex, date_diff, obs_diff, age) %>%
group_by(ring_ID) %>%
dplyr::summarise(mean_obs_diff = mean(obs_diff, na.rm = TRUE),
median_obs_diff = median(obs_diff, na.rm = TRUE)) %>%
arrange(desc(median_obs_diff)) %>%
filter(median_obs_diff <= 2)
# subset to only include individuals with mainly observations spaced 2 days apart
coucal_wp_clean <-
coucal_wp_clean %>%
filter(ring_ID %in% two_day_rings$ring_ID) %>%
# remove questionable observation 18 km
filter(distance < 18000)#### Modeling sex-specific movements ####
# Model 1 ####
# Does daily distance moved increase with age since leaving the nest
# and does this relationship differ between males and females?
# M1 Black Coucals: strong effect of age, no interaction with sex
mod_BC_dist <-
lmer(log(distance) ~ poly(age, 2) * sex + (1|ring_ID),
data = filter(coucal_wp_clean, species == "BC"))
# effect-size plots
coefplot::coefplot(mod_BC_dist)
plot(allEffects(mod_BC_dist))
# extract fitted values
mod_BC_dist_fits <-
as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_dist,
xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
sex = c("Female", "Male"))))
mod_BC_dist_fits$species <- "BC"
# M1 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_dist <-
lmer(log(distance) ~ poly(age, 2) * sex +
(1|ring_ID),
data = filter(coucal_wp_clean, species == "WBC" & distance > 0))
# effect-size plots
coefplot::coefplot(mod_WBC_dist)
plot(allEffects(mod_WBC_dist))
# extract fitted values
mod_WBC_dist_fits <-
as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_dist,
xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
sex = c("Female", "Male"))))
mod_WBC_dist_fits$species <- "WBC"
mod_dist_fits <-
bind_rows(mod_BC_dist_fits, mod_WBC_dist_fits)# first set up plotting parameters
# define the plotting theme to be used in subsequent ggplots
luke_theme <-
theme_bw() +
theme(
text = element_text(family = "Verdana"),
legend.title = element_text(size = 16),
legend.text = element_text(size = 12),
axis.title.x = element_text(size = 12),
axis.text.x = element_text(size = 12),
axis.title.y = element_text(size = 16),
axis.text.y = element_text(size = 12),
strip.text = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size = 0.5, colour = "grey40"),
axis.ticks.length = unit(0.2, "cm"),
panel.border = element_rect(linetype = "solid", colour = "grey"),
legend.position = c(0.1, 0.9)
)
# set plotting color palettes
plot_palette_sex <- RColorBrewer::brewer.pal(8, "Dark2")[c(2,1)]
plot_palette_brood <- RColorBrewer::brewer.pal(7, "Blues")[c(3:7)]
# specify the facet labels for each species
species.labs <- c("Black Coucal", "White-browed Coucal")
names(species.labs) <- c("BC", "WBC")# Model 2 ####
# Does cumulative distance moved increase with age since leaving the nest
# and does this relationship differ between males and females?
# M2 Black Coucals: strong effect of age, no interaction with sex
mod_BC_cum_dist <-
lmer(log(cum_distance) ~ poly(age, 2) * sex + (1|ring_ID),
data = filter(coucal_wp_clean, species == "BC"))
# effect-size plots
coefplot::coefplot(mod_BC_cum_dist)
plot(allEffects(mod_BC_cum_dist))
# extract fitted values
mod_BC_cum_dist_fits <-
as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_cum_dist,
xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
sex = c("Female", "Male"))))
mod_BC_cum_dist_fits$species <- "BC"
# M2 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_cum_dist <-
lmer(log(cum_distance) ~ poly(age, 2) * sex +
(1|ring_ID),
data = filter(coucal_wp_clean, species == "WBC" & distance > 0))
# effect-size plots
coefplot::coefplot(mod_WBC_cum_dist)
plot(allEffects(mod_WBC_cum_dist))
# extract fitted values
mod_WBC_cum_dist_fits <-
as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_cum_dist,
xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
sex = c("Female", "Male"))))
mod_WBC_cum_dist_fits$species <- "WBC"
mod_cum_dist_fits <-
bind_rows(mod_BC_cum_dist_fits, mod_WBC_cum_dist_fits)# Model 3 ####
# Do individuals move further from the nest as they age
# and does this relationship differ between males and females?
# M3 Black Coucals: strong effect of age, no interaction with sex
mod_BC_disp <-
lmer(log(distance_origin) ~ poly(age, 2) * sex + (1|ring_ID),
data = filter(coucal_wp_clean, species == "BC" & distance_origin > 0))
# effect-size plots
coefplot::coefplot(mod_BC_disp)
plot(allEffects(mod_BC_disp))
# extract fitted values
mod_BC_disp_fits <-
as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_BC_disp,
xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]),
max(coucal_wp_clean[coucal_wp_clean$species == "BC", "age"]), 1),
sex = c("Female", "Male"))))
mod_BC_disp_fits$species <- "BC"
# M3 White-browed Coucals: strong effect of age, no interaction with sex
mod_WBC_disp <-
lmer(log(distance_origin) ~ poly(age, 2) * sex +
(1|ring_ID),
data = filter(coucal_wp_clean, species == "WBC" & distance_origin > 0))
# effect-size plots
coefplot::coefplot(mod_WBC_disp)
plot(allEffects(mod_WBC_disp))
# extract fitted values
mod_WBC_disp_fits <-
as.data.frame(effect(term = "poly(age, 2) * sex", mod = mod_WBC_disp,
xlevels = list(age = seq(min(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]),
max(coucal_wp_clean[coucal_wp_clean$species == "WBC", "age"]), 1),
sex = c("Female", "Male"))))
mod_WBC_disp_fits$species <- "WBC"
mod_disp_fits <-
bind_rows(mod_BC_disp_fits, mod_WBC_disp_fits)# M1 plot
ggplot() +
luke_theme +
geom_line(data = coucal_wp_clean,
aes(x = age, y = log(distance), group = ring_ID, color = sex),
alpha = 0.2) +
geom_line(data = mod_dist_fits, aes(x = age, y = fit, color = sex),
lwd = 0.5) +
geom_ribbon(data = mod_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
lwd = 1, alpha = 0.25) +
facet_grid(species ~ .,
labeller = labeller(species = species.labs)) +
ylab(expression(paste("Distance moved between observations (log10)" %+-% "95% CI", sep = ""))) +
xlab("Age since leaving nest") +
scale_color_manual(values = plot_palette_sex,
name = "Sex",
breaks = c("Female", "Male"),
labels = c("Female", "Male")) +
scale_fill_manual(values = plot_palette_sex,
name = "Sex",
breaks = c("Female", "Male"),
labels = c("Female", "Male"))
# M2 plot
ggplot() +
luke_theme +
geom_line(data = coucal_wp_clean,
aes(x = age, y = log(cum_distance), group = ring_ID, color = sex),
alpha = 0.2) +
geom_line(data = mod_cum_dist_fits, aes(x = age, y = fit, color = sex),
lwd = 0.5) +
geom_ribbon(data = mod_cum_dist_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
lwd = 1, alpha = 0.25) +
facet_grid(species ~ .,
labeller = labeller(species = species.labs)) +
ylab(expression(paste("Cumulative distance moved (log10)" %+-% "95% CI", sep = ""))) +
xlab("Age since leaving nest") +
scale_color_manual(values = plot_palette_sex,
name = "Sex",
breaks = c("Female", "Male"),
labels = c("Female", "Male")) +
scale_fill_manual(values = plot_palette_sex,
name = "Sex",
breaks = c("Female", "Male"),
labels = c("Female", "Male"))
# M3 plot
ggplot() +
luke_theme +
geom_line(data = filter(coucal_wp_clean, distance_origin > 0),
aes(x = age, y = log(distance_origin), group = ring_ID, color = sex),
alpha = 0.2) +
geom_line(data = mod_disp_fits, aes(x = age, y = fit, color = sex),
lwd = 0.5) +
geom_ribbon(data = mod_disp_fits, aes(x = age, ymax = upper, ymin = lower, fill = sex),
lwd = 1, alpha = 0.25) +
facet_grid(species ~ .,
labeller = labeller(species = species.labs)) +
ylab(expression(paste("Daily distance from nest (log10)" %+-% "95% CI", sep = ""))) +
xlab("Age since leaving nest") +
scale_color_manual(values = plot_palette_sex,
name = "Sex",
breaks = c("Female", "Male"),
labels = c("Female", "Male")) +
scale_fill_manual(values = plot_palette_sex,
name = "Sex",
breaks = c("Female", "Male"),
labels = c("Female", "Male"))
Supplementary Analysis 2: weekly sex-specific adult survival
status_dat_ad <-
# read raw data
read.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv",
header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
# rename ring_ID column
dplyr::rename(ring_ID = Alu,
ageC2 = days_from_ringing_until_last_observation) %>%
# specify date strings properly
mutate(month_ringed = str_sub(date_dec_initially_ringed, start = 5, end = 6),
day_ringed = str_sub(date_dec_initially_ringed, start = 7, end = 8),
year_ringed = str_sub(date_dec_initially_ringed, start = 1, end = 4),
month_last_ob = str_sub(date_dec_last_observed, start = 5, end = 6),
day_last_ob = str_sub(date_dec_last_observed, start = 7, end = 8),
year_last_ob = str_sub(date_dec_initially_ringed, start = 1, end = 4)) %>%
mutate(date_ringed = as.Date(paste(year_ringed, month_ringed, day_ringed, sep = "-"),
format = "%Y-%m-%d"),
date_last_ob = as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep = "-"),
format = "%Y-%m-%d")) %>%
# select variables of interest
dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>%
# remove all white space from data
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
# specify empty data as NA
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
# exclude all individuals that died in the nest
# filter(Fledged_status == "yes") %>%
# classify columns
mutate(sex = as.factor(sex),
ageC2 = as.numeric(ageC2)) %>%
# remove rows with missing sex, age, and status info
filter(!is.na(sex) & !is.na(ageC2) & !is.na(statusC)) %>%
# make a unique id for each individual
mutate(# create the age of entry into the data (all at age 15)
entry = 0,
# specify the age of death or censoring
exit = ageC2,
# make the event numeric and specify if
# the individual died (1) or was censored (0)
event = as.numeric(statusC),
sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex))) %>%
# consolidate to variables of interest
dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob)
status_dat_ad_week <-
status_dat_ad %>%
dplyr::select(species, ring_ID, sex, date_ringed, date_last_ob) %>%
pivot_longer(-c(species, ring_ID, sex), names_to = "status", values_to = "date") %>%
mutate(year_numeric = as.numeric(lubridate::year(date)),
week_numeric = as.numeric(strftime(date, format = "%V"))) %>%
dplyr::select(species, ring_ID, year_numeric, sex, week_numeric)
# import raw csv data into R
BC_detect_dat_A <-
read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>%
dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>%
mutate(month = str_sub(date_dec, start = 5, end = 6),
day = str_sub(date_dec, start = 7, end = 8)) %>%
mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>%
mutate(across(c("sex", "site", "age_status"), tolower)) %>%
mutate(age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>%
mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>%
mutate(sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex)),
month_year = format(date, "%Y-%m"),
month_numeric = as.numeric(month),
year_numeric = as.numeric(year),
week_numeric = as.numeric(strftime(date, format = "%V"))) %>%
filter(species == "BC") %>%
filter(age_status == "A") %>%
dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%
bind_rows(., filter(status_dat_ad_week, species == "BC")) %>%
filter(!is.na(week_numeric)) %>%
distinct() %>%
mutate(week_std = round(scale_by(week_numeric ~ year_numeric, ., scale = 0), digits = 0)) %>%
mutate(week_std = week_std + abs(min(week_std, na.rm = TRUE))) %>%
dplyr::rename(year = year_numeric)
WBC_detect_dat_A <-
read_xlsx("data/raw/All_coucal_waypoints_2001_2019_20200202.xlsx", na = "NA", col_types = "text") %>%
dplyr::select(species, ring_ID, sex, year, site, age_status, date_dec) %>%
mutate(month = str_sub(date_dec, start = 5, end = 6),
day = str_sub(date_dec, start = 7, end = 8)) %>%
mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>%
mutate(across(c("sex", "site", "age_status"), tolower)) %>%
mutate(age_status = ifelse(age_status == "adult", "A", ifelse(age_status == "juvenile", "J", "XXX")),
ring_ID = str_replace_all(string = ring_ID, fixed(" "), "")) %>%
mutate(across(c("species", "ring_ID", "sex", "site", "age_status"), as.factor)) %>%
mutate(sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex)),
month_year = format(date, "%Y-%m"),
month_numeric = as.numeric(month),
year_numeric = as.numeric(year),
week_numeric = as.numeric(strftime(date, format = "%V"))) %>%
filter(species == "WBC") %>%
filter(age_status == "A") %>%
dplyr::select(species, ring_ID, year_numeric, sex, week_numeric) %>%
bind_rows(., filter(status_dat_ad_week, species == "WBC")) %>%
filter(!is.na(week_numeric)) %>%
distinct() %>%
mutate(week_std = round(scale_by(week_numeric ~ year_numeric, ., scale = 0), digits = 0)) %>%
mutate(week_std = week_std + abs(min(week_std, na.rm = TRUE))) %>%
dplyr::rename(year = year_numeric)
# use the BaSTA function "CensusToCaptHist()" function to convert long format
# encounter histories of each individual, to wide format with 1's and 0's for
# each week of encounter for each year
BC_weekly_A_ch <-
CensusToCaptHist(ID = BC_detect_dat_A$ring_ID,
d = BC_detect_dat_A$week_std,
dformat = "%W", timeInt = "%W") %>%
dplyr::rename(ring_ID = ID) %>%
mutate(ID = rownames(.)) %>%
left_join(., dplyr::select(BC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>%
distinct()
WBC_weekly_A_ch <-
CensusToCaptHist(ID = WBC_detect_dat_A$ring_ID,
d = WBC_detect_dat_A$week_std,
dformat = "%W", timeInt = "%W") %>%
dplyr::rename(ring_ID = ID) %>%
mutate(ID = rownames(.)) %>%
left_join(., dplyr::select(WBC_detect_dat_A, ring_ID, sex, year), by = "ring_ID") %>%
distinct()
# create a two-character string for each encounter and clean the output
BC_weekly_A_ch_Burnham <-
sapply(BC_weekly_A_ch[2:26], function(x) paste(x, "0", sep = "")) %>%
cbind(BC_weekly_A_ch[c(1,45:length(BC_weekly_A_ch))]) %>%
mutate(across(everything(), ~str_replace(string = .x, pattern = "NA0", "00"))) %>%
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
distinct() %>%
mutate(sex = as.factor(sex))
WBC_weekly_A_ch_Burnham <-
sapply(WBC_weekly_A_ch[2:26], function(x) paste(x, "0", sep = "")) %>%
cbind(WBC_weekly_A_ch[c(1,54:length(WBC_weekly_A_ch))]) %>%
mutate(across(everything(), ~str_replace(string = .x, pattern = "NA0", "00"))) %>%
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
distinct() %>%
mutate(sex = as.factor(sex))
# Import status data
status_dat_all_ad <-
# read raw data
read.csv("data/raw/Coucal_adults_survival_2001-2019_20200129.csv",
header = TRUE, stringsAsFactors = FALSE, na.strings = c("", " ", "NA")) %>%
# rename ring_ID column
dplyr::rename(ring_ID = Alu,
ageC2 = days_from_ringing_until_last_observation) %>%
# specify date strings properly
mutate(month_ringed = str_sub(date_dec_initially_ringed, start = 5, end = 6),
day_ringed = str_sub(date_dec_initially_ringed, start = 7, end = 8),
year_ringed = str_sub(date_dec_initially_ringed, start = 1, end = 4),
month_last_ob = str_sub(date_dec_last_observed, start = 5, end = 6),
day_last_ob = str_sub(date_dec_last_observed, start = 7, end = 8),
year_last_ob = str_sub(date_dec_initially_ringed, start = 1, end = 4)) %>%
mutate(date_ringed = as.Date(paste(year_ringed, month_ringed, day_ringed, sep = "-"),
format = "%Y-%m-%d"),
date_last_ob = as.Date(paste(year_last_ob, month_last_ob, day_last_ob, sep = "-"),
format = "%Y-%m-%d")) %>%
# select variables of interest
dplyr::select(species, ring_ID, sex, year, ageC2, ageC, statusC, date_ringed, date_last_ob) %>%
# remove all white space from data
mutate(across(everything(), ~str_trim(.x))) %>%
mutate(across(everything(), ~str_replace_all(.x, fixed(" "), ""))) %>%
# specify empty data as NA
mutate(across(everything(), ~gsub("^$|^ $", NA, .x))) %>%
# classify columns
mutate(sex = as.factor(sex),
ageC2 = as.numeric(ageC2)) %>%
# remove rows with missing sex, age, and status info
filter(!is.na(sex) & !is.na(ageC2) & !is.na(statusC)) %>%
# make a unique id for each individual
mutate(# create the age of entry into the data (all at age 15)
entry = 0,
# specify the age of death or censoring
exit = ageC2,
# make the event numeric and specify if
# the individual died (1) or was censored (0)
event = as.numeric(statusC),
sex = ifelse(str_detect(sex, pattern = "[Ff]emale"), "F",
ifelse(str_detect(sex, pattern = "[Mm]ale"), "M", sex))) %>%
# consolidate to variables of interest
dplyr::select(species, ring_ID, year, sex, entry, exit, event, date_ringed, date_last_ob) %>%
dplyr::select(species, ring_ID, year, sex, event)
# join detection history with status data
BC_detect_status_join <-
filter(status_dat_all_ad, species == "BC") %>%
left_join(BC_weekly_A_ch_Burnham, ., by = c("ring_ID", "year", "sex")) %>%
dplyr::select(-species) %>%
mutate(event = ifelse(is.na(event), 0, event))
WBC_detect_status_join <-
filter(status_dat_all_ad, species == "WBC") %>%
left_join(WBC_weekly_A_ch_Burnham, ., by = c("ring_ID", "year", "sex")) %>%
dplyr::select(-species) %>%
mutate(event = ifelse(is.na(event), 0, event))
# determine the last and first detection for each individual
max_age_index_BC <-
apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>%
lapply(., function(x) x[which.max(x)]) %>%
unlist(.)
min_age_index_BC <-
apply(BC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>%
lapply(., function(x) x[which.min(x)]) %>%
unlist(.)
max_age_index_WBC <-
apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>%
lapply(., function(x) x[which.max(x)]) %>%
unlist(.)
min_age_index_WBC <-
apply(WBC_detect_status_join[, c(1:25)], 1, function(x) which(x == "10")) %>%
lapply(., function(x) x[which.min(x)]) %>%
unlist(.)
# put "11" at the last detection if the status was a "1" (i.e., dead)
for(i in 1:nrow(BC_detect_status_join)){
if(BC_detect_status_join$event[i] == "1"){
BC_detect_status_join[i, max_age_index_BC[i]] <- "11"
}
}
for(i in 1:nrow(WBC_detect_status_join)){
if(WBC_detect_status_join$event[i] == "1"){
WBC_detect_status_join[i, max_age_index_WBC[i]] <- "11"
}
}
# bind the min and max age info to the capture history data
detect_status_BC <-
bind_cols(BC_detect_status_join, max_age_index_BC) %>%
bind_cols(., min_age_index_BC) %>%
dplyr::rename(std_week_last_obs = '...31',
std_week_first_obs = '...32')
detect_status_WBC <-
bind_cols(WBC_detect_status_join, max_age_index_WBC) %>%
bind_cols(., min_age_index_WBC) %>%
dplyr::rename(std_week_last_obs = '...31',
std_week_first_obs = '...32')
# consolidate capture histories for each sex
Black_Coucal_adult_weekly_Burnham_ch <-
apply(detect_status_BC[, c(1:25)], 1, paste, collapse = "") %>%
as.data.frame() %>%
dplyr::rename(ch = '.') %>%
mutate(ch = as.character(ch)) %>%
dplyr::bind_cols(., detect_status_BC[, c(26:length(detect_status_BC))]) %>%
filter(str_detect(ch, "11", negate = TRUE) | str_count(ch, "1") > 2) %>%
dplyr::select(ch, ring_ID, sex, year) %>%
drop_na() %>%
mutate(year = ifelse(year == "1006", "2006", year)) %>%
mutate_at(vars(ring_ID, sex, year), as.factor)
White_browed_Coucal_adult_weekly_Burnham_ch <-
apply(detect_status_WBC[, c(1:25)], 1, paste, collapse = "") %>%
as.data.frame() %>%
dplyr::rename(ch = '.') %>%
mutate(ch = as.character(ch)) %>%
dplyr::bind_cols(., detect_status_WBC[, c(26:length(detect_status_WBC))]) %>%
filter(str_detect(ch, "11", negate = TRUE) | str_count(ch, "1") > 2) %>%
dplyr::select(ch, ring_ID, sex, year) %>%
drop_na() %>%
mutate_at(vars(ring_ID, sex, year), as.factor)
BC_adult_ch <-
Black_Coucal_adult_weekly_Burnham_ch %>%
mutate(sex = as.factor(sex),
year = as.factor(year)) %>%
dplyr::select(ch, sex, year)
WBC_adult_ch <-
White_browed_Coucal_adult_weekly_Burnham_ch %>%
mutate(sex = as.factor(sex),
year = as.factor(year)) %>%
dplyr::select(ch, sex, year)
BC_adult_ch %>%
mutate(status = ifelse(str_detect(ch, "11"), "dead", "censored")) %>%
group_by(sex, status) %>%
dplyr::summarise(n_ind = n())
WBC_adult_ch %>%
mutate(status = ifelse(str_detect(ch, "11"), "dead", "censored")) %>%
group_by(sex, status) %>%
dplyr::summarise(n_ind = n())
length(levels(BC_adult_ch$year))
length(levels(WBC_adult_ch$year))# process the adult data as a "Burnham" analysis
BC_coucal_adult.proc <- RMark::process.data(data = BC_adult_ch,
model = "Burnham",
groups = c("sex", "year"))
WBC_coucal_adult.proc <- RMark::process.data(data = WBC_adult_ch,
model = "Burnham",
groups = c("sex", "year"))
# make design matrix
BC_coucal_adult.ddl <-
RMark::make.design.data(BC_coucal_adult.proc)
WBC_coucal_adult.ddl <-
RMark::make.design.data(WBC_coucal_adult.proc)
burnham_S_analysis_run_pFr_dot <-
function(proc_data, design_data){
# sex-specific model
S.sex = list(formula = ~sex)
# sex- and week-specific model
S.sex_Time = list(formula = ~sex + Time)
# sex- and week-specific model (quadratic)
S.sex_Quad = list(formula = ~sex + (I(Time) + I(Time^2)))
# sex- and week-specific model (interaction)
S.sex_x_Time = list(formula = ~sex * Time)
# sex- and week-specific model (quadratic with interaction)
S.sex_x_Quad = list(formula = ~sex * (I(Time) + I(Time^2)))
# p (encounter probability):
# constant model
p.dot = list(formula = ~1)
# F (site fidelity probability):
# constant model
F.dot = list(formula = ~1)
# r (recovery probability)
# constant model
r.dot = list(formula = ~1)
# create model list for all a priori models above
cml <- RMark::create.model.list("Burnham")
# run the models in program MARK
model.list <- RMark::mark.wrapper(cml, data = proc_data,
ddl = design_data, delete = TRUE,
wrap = FALSE, threads = 1, brief = TRUE,
silent = TRUE, output = FALSE)
# output the model list and sotre the results
return(model.list)
}
adult_S_analysis_out_BC <-
burnham_S_analysis_run_pFr_dot(proc_data = BC_coucal_adult.proc,
design_data = BC_coucal_adult.ddl)
adult_S_analysis_out_WBC <-
burnham_S_analysis_run_pFr_dot(proc_data = WBC_coucal_adult.proc,
design_data = WBC_coucal_adult.ddl)
BC_sex_x_quad_reals <-
adult_S_analysis_out_BC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>%
bind_cols(data.frame(str_split_fixed(rownames(.), " ",
n = 5)), .) %>%
mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F",
"Female", "Male"))) %>%
filter(X1 == "S") %>%
dplyr::select(sex, age, estimate, se, lcl, ucl) %>%
mutate(species = "Black Coucal")
WBC_sex_x_quad_reals <-
adult_S_analysis_out_WBC$S.sex_x_Quad.p.dot.r.dot.F.dot$results$real %>%
bind_cols(data.frame(str_split_fixed(rownames(.), " ",
n = 5)), .) %>%
mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F",
"Female", "Male"))) %>%
filter(X1 == "S") %>%
dplyr::select(sex, age, estimate, se, lcl, ucl) %>%
mutate(species = "White-browed Coucal")
BC_sex_reals <-
adult_S_analysis_out_BC$S.sex.p.dot.r.dot.F.dot$results$real %>%
bind_cols(data.frame(str_split_fixed(rownames(.), " ",
n = 5)), .) %>%
mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F",
"Female", "Male"))) %>%
filter(X1 == "S") %>%
dplyr::select(sex, age, estimate, se, lcl, ucl) %>%
mutate(species = "Black Coucal")
WBC_sex_reals <-
adult_S_analysis_out_WBC$S.sex.p.dot.r.dot.F.dot$results$real %>%
bind_cols(data.frame(str_split_fixed(rownames(.), " ",
n = 5)), .) %>%
mutate(age = as.integer(unlist(str_extract_all(X4,"[0-9]+"))),
sex = as.factor(ifelse(unlist(str_extract_all(X2,"[FM]")) == "F",
"Female", "Male"))) %>%
filter(X1 == "S") %>%
dplyr::select(sex, age, estimate, se, lcl, ucl) %>%
mutate(species = "White-browed Coucal")
time_varying_plot <-
ggplot(data = bind_rows(BC_sex_x_quad_reals, WBC_sex_x_quad_reals)) +
geom_line(aes(x = age, y = estimate, group = sex, color = sex), lwd = 0.8) +
geom_ribbon(aes(x = age, ymin = lcl, ymax = ucl, group = sex, fill = sex), alpha = 0.5) +
scale_color_manual(values = sex_pal2,
labels = c("Female", "Male")) +
scale_fill_manual(values = sex_pal2,
labels = c("Female", "Male")) +
facet_grid(. ~ species) +
theme_bw() +
theme(text = element_text(family = "Verdana"),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
strip.text = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size = 0.5, colour = "grey40"),
axis.ticks.length = unit(0.2, "cm"),
panel.border = element_rect(linetype = "solid", colour = "grey"),
legend.position = c(0.8, 0.2),
legend.title = element_blank()) +
xlab("Week of breeding season (standardized across years)") +
ylab("Adult weekly survival probability (± 95% CI)")
constant_plot <-
ggplot2::ggplot() +
geom_errorbar(data = bind_rows(BC_sex_reals, WBC_sex_reals),
aes(x = sex, ymax = ucl, ymin = lcl, color = sex),
alpha = 1, width = 0.05, lwd = 0.5) +
geom_point(data = bind_rows(BC_sex_reals, WBC_sex_reals),
aes(x = sex, y = estimate, fill = sex),
lwd = 3, shape = 21) +
scale_color_manual(values = sex_pal2,
labels = c("Female", "Male")) +
scale_fill_manual(values = sex_pal2,
labels = c("Female", "Male")) +
facet_grid(. ~ species) +
theme_bw() +
theme(
text = element_text(family = "Verdana"),
axis.title = element_blank(),
axis.text = element_text(size = 10),
strip.text = element_text(size = 12),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks = element_line(size = 0.5, colour = "grey40"),
axis.ticks.length = unit(0.2, "cm"),
panel.border = element_rect(linetype = "solid", colour = "grey"),
legend.position = "none") +
ylab("Adult weekly survival probability (± 95% CI)") +
scale_y_continuous(limits = c(0, 1))
combo_plot_S <-
time_varying_plot + constant_plot + plot_annotation(tag_levels = "A")
# ggsave(combo_plot_S,
# filename = "products/figures/combo_plot_S_weekly.jpeg",
# width = 12,
# height = 6,
# units = "in",
# dpi = 300)Supplementary Analysis 3: Plumage-based sex-specific age structure
age_plumage_dist_data <-
read_xlsx("data/raw/age_classes_2001-2020_2.xlsx", na = "NA", col_types = "text") %>%
dplyr::rename_all(~str_replace_all(., "\\s+", "")) %>%
dplyr::mutate(species = tolower(species)) %>%
dplyr::filter(str_detect(string = species, pattern = "black")) %>%
dplyr::mutate(month = str_sub(`date(yymmdd)`, start = 5, end = 6),
day = str_sub(`date(yymmdd)`, start = 7, end = 8)) %>%
dplyr::mutate(date = as.Date(paste(year, month, day, sep = "-"), format = "%Y-%m-%d")) %>%
# remove all white space from data
dplyr::mutate(dplyr::across(.cols = dplyr::everything(),
str_remove_all, pattern = fixed(" "))) %>%
dplyr::mutate(species = "BC",
sex = ifelse(sex == "female", "F", ifelse(sex == "male", "M", "XXX"))) %>%
dplyr::mutate(dplyr::across(.cols = dplyr::everything(),
str_replace_all, pattern = fixed("onebarred"), replacement = fixed("barred"))) %>%
dplyr::mutate(Secondarycoverts = ifelse(SN == "462", "rufous", Secondarycoverts)) %>%
dplyr::mutate(Secondarycoverts = ifelse(SN == "409", "barred", Secondarycoverts)) %>%
dplyr::mutate(age_LEP = ifelse(str_detect(primarycoverts, "barred") &
str_detect(primaries, "barred") &
str_detect(Secondarycoverts, "barred") &
str_detect(secondaries, "barred"), 1,
ifelse(str_detect(primaries, "barred") |
str_detect(secondaries, "barred"), 2,
ifelse((str_detect(primarycoverts, "barred") |
str_detect(Secondarycoverts, "barred")) &
(str_detect(primaries, "rufous") |
str_detect(secondaries, "rufous")), 2,
ifelse(str_detect(primarycoverts, "rufous") &
str_detect(primaries, "rufous") &
str_detect(Secondarycoverts, "rufous") &
str_detect(secondaries, "rufous"), 3, 0))))) %>%
dplyr::mutate(age = ifelse(age == "juvenile", 1, age)) %>%
dplyr::mutate(CHECK = ifelse(age != age_LEP, "CHECK", "")) %>%
dplyr::mutate(age_LEP = age) %>%
dplyr::filter(!is.na(age_LEP)) %>%
dplyr::filter(year != "2020") %>%
dplyr::group_by(year, sex, age_LEP) %>%
dplyr::summarise(n_ages = n_distinct(Alu)) %>%
dplyr::group_by(year, sex) %>%
dplyr::mutate(prop = n_ages/sum(n_ages)) %>%
dplyr::mutate(age_LEP = as.character(age_LEP)) %>%
dplyr::mutate(age_LEP = ifelse(age_LEP == "3", "3+", age_LEP)) %>%
dplyr::mutate(age_LEP = fct_rev(as.factor(age_LEP)))
age_plumage_dist_data2 <-
age_plumage_dist_data %>%
dplyr::group_by(sex, age_LEP) %>%
dplyr::summarise(mean_n_ages = mean(n_ages)) %>%
dplyr::group_by(sex) %>%
dplyr::mutate(prop = mean_n_ages/sum(mean_n_ages),
year = "Grand\naverage") %>%
dplyr::bind_rows(age_plumage_dist_data, .) %>%
dplyr::mutate(sex = ifelse(sex == "F", "Female", "Male"))
age_structure_plot <-
ggplot(age_plumage_dist_data2,
aes(fill = sex, alpha = age_LEP,
y = prop, x = sex)) +
geom_bar(position="stack", stat="identity") +
facet_wrap(. ~ year) +
theme_bw() +
theme(text = element_text(family = "Verdana"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.ticks.length = unit(0.1, "cm"),
panel.border = element_blank(),
panel.spacing = unit(0.3, "lines"),
strip.background = element_blank(),
legend.position = "top") +
scale_fill_manual(values = sex_pal2, name = "Sex", guide = FALSE) +
scale_alpha_manual(values = c(0.3, 0.65, 1), name = "Age (years)",
guide = guide_legend(title.position = "top", title.hjust = 0.5, label.position = "bottom")) +
ylab("Proportion of breeding adults") +
xlab("Sex")
# ggsave(age_structure_plot,
# filename = "products/figures/age_structure_plot.jpeg",
# width = 6,
# height = 6,
# units = "in",
# dpi = 300)
age_structure_table <-
age_plumage_dist_data2 %>%
dplyr::filter(str_detect(year, "Grand")) %>%
dplyr::ungroup() %>%
dplyr::select(sex, age_LEP, mean_n_ages, prop) %>%
dplyr::arrange(sex, desc(age_LEP)) %>%
gt(rowname_col = "row",
groupname_col = "sex") %>%
cols_label(age_LEP = html("<i>Age class</i>"),
mean_n_ages = "Annual\nmean no. ind.",
prop = "Proportion") %>%
fmt_number(columns = vars(mean_n_ages, prop),
decimals = 2,
use_seps = FALSE) %>%
cols_align(align = "right",
columns = vars(mean_n_ages)) %>%
cols_align(align = "left",
columns = vars(prop)) %>%
tab_options(row_group.font.weight = "bold",
row_group.background.color = brewer.pal(9,"Greys")[3],
table.font.size = 12,
data_row.padding = 3,
row_group.padding = 4,
summary_row.padding = 2,
column_labels.font.size = 14,
row_group.font.size = 12,
table.width = pct(60))
# ggsave(age_structure_table,
# filename = "products/tables/age_structure_table.jpeg",
# width = 4,
# height = 4,
# units = "in",
# dpi = 300)